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Abstract. We introduce a new approach to high-order accuracy for the numerical solution of diffusion 
problems by solving the equations in differential form using a reconstruction technique. The approach has 
the advantages of simplicity and economy. It results in several new high-order methods including a 
simplified version of discontinuous Galerkin (DG). It also leads to new definitions of common value and 
common gradient — quantities at each interface shared by the two adjacent cells. In addition, the new 
approach clarifies the relations among the various choices of new and existing common quantities. Fourier 
stability and accuracy analyses are carried out for the resulting schemes. Extensions to the case of 
quadrilateral meshes are obtained via tensor products. For the two-point boundary value problem (steady 
state), it is shown that these schemes, which include most popular DG methods, yield exact common 
interface quantities as well as exact cell average solutions for nearly all cases. 
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1. Introduction 

In the field of Computational Fluid Dynamics, low-order methods are less accurate, but generally are 
robust and reliable; therefore, they are routinely employed in practical calculations. High-order (third and 
above) methods can provide accurate solutions; however, they are more complicated and less robust. The 
need to improve and simplify high-order methods as well as to develop new ones with more favorable 
properties has attracted the interest of many researchers. The current work is a step in this direction. 

For simplicity, we discuss the case of one spatial dimension, which contains all of the key new ideas. 
For this case, the solution u is approximated in each cell (interval) by the data at K points called solution 
points ( K > 1 ). These K pieces of data determine a polynomial of degree K— 1 , herein called a solution 
polynomial. Such polynomials collectively form a function, which is generally discontinuous across cell 
interfaces. One advantage of allowing such discontinuities is that the resulting method is local. Another 
advantage is that the method can better resolve discontinuities such as shocks in the exact solution. 
Calculating the first and second derivatives of u by those of the solution polynomials leads to incorrect 
results, since these quantities involve no interaction of the data between cells. For interaction to take place, 
at each interface, we need to define a value u and a derivative (gradient) value u x common to the two 
adjacent cells; the first and second derivative estimates making use of these common quantities (value or 
derivative) would involve interaction. These common quantities can be employed in two different 
approaches. 

The first approach, introduced by Reed and Hill (1973) for neutron transport equations, is the 
discontinuous Galerkin or DG approach. It was developed for fluid dynamics equations by Cockburn and 
Shu (1989), Bassi and Rebay (1997a, b), and numerous other researchers; see, e.g., the review paper by 
Cockburn, Karniadakis, and Shu (2000). Aside from the discontinuous representation of the solution, 
another key ingredient of the DG methods is the integral formulation using test functions. The DG 
approach was applied to solve diffusion-type problems by several authors, e.g., Arnold (1982), Bassi and 
Rebay (1997a), Cockburn and Shu (1998), Oden et al. (1998), Arnold et al. (2002). The common derivative 
at each interface, which leads to a scheme with a rather wide stencil called BR1, was presented in (Bassi 
and Rebay 1997a). It was modified so that the resulting scheme has a compact stencil including only 
immediate neighbors (Bassi and Rebay 2000). This modification, often called BR2, was analyzed in 
(Brezzi et al. 2000). In addition to the advantage of having a compact stencil, which simplifies the coding 
of boundary conditions, another advantage of the BR2 scheme is that for the steady state case, the 
coefficient matrix for the solution is invertible. The technique of compact stencil was applied by Peraire 
and Persson (2008) to the case of one-sided common quantities; the resulting method, which is a 
modification of the local DG or LDG scheme of Cockburn and Shu (1998), is named compact DG or CDG. 
Finally, Van Leer and Nomura (2006) obtained common quantities by a recovery principle: a polynomial of 
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degree 2 K - 1 approximating the solution polynomials on a domain formed by the two cells adjacent to the 
interface is constructed via least squares. 

The second approach, introduced in (Kopriva and Kolias 1996) and (Kopriva 1998) for a quadrilateral 
mesh, solves the differential form of the equation as opposed to the integral form. Here, the solutions are 
evaluated at the K Chebyshev-Gauss points, and the fluxes, at the K + 1 Chebyshev-Lobatto points. For 
diffusion problems, the common quantities are defined as in the BR1 scheme. This approach on a triangular 
mesh results in a scheme called spectral difference (Liu et al. 2006). These solution-point and flux-point (or 
staggered-mesh) methods have the advantage of conceptual simplicity and the disadvantages of using two 
sets of points (or grids) as well as being mildly unstable (Huynh 2007, Van den Abeele et al. 2008). 

In this paper, we present a new approach to high-order accuracy for diffusion problems using the 
differential form. The approach has the advantage of simplicity and economy. It employs one set of points, 
namely, the solution points, and fluxes are evaluated only at the cell boundaries. It was introduced to solve 
advection-type problems in (Huynh 2007), where the objective was the evaluation of the first derivative. 
Here, the approach is applied to solve diffusion problems, and the objective is the evaluation of the second 
derivative. To this end, first, the common value at each interface can be defined as the average of the two 
values just to its left and right as in the BR1 and BR2 schemes, or a one-sided value, say, the one to its left 
as in the LDG and CDG scheme (other criteria will be introduced later). Next, we reconstruct the solution 
by a function which is continuous across the interfaces and, in each cell, is a polynomial of one degree 
higher than that of the solution polynomial so that its derivative is of the same degree. To assure continuity 
across the interfaces, this reconstruction is required to take on the common interface value. In each cell, the 
reconstruction is obtained by adding two correction functions to the solution polynomial, one to match the 
common value at the left interface, and the other to match the common value at the right. Due to symmetry, 
we only need to focus on the correction for the left interface. By rescaling, the problem reduces to defining 
a correction function g on the interval [-1,1] such that g(-l) = l, g(l) = 0, and g is a polynomial of 

degree K approximating the zero function in some sense. The condition g(-l) = 1 deals with the jump at 
the left interface whereas the condition g(l) = 0 leaves the right interface value unchanged. The crucial 
condition of approximating the zero function can be met by the least-squares principle: g is required to be 
orthogonal to all polynomials of degree K — 2 or less. Such a correction function turns out to be the Radau 
polynomial, and the resulting method is identical to DG. The current version of the scheme, however, is 
simpler than the standard version using integral formulation: at each solution point, it involves only a 
straightforward derivative formula and a simple correction term regardless of the choice of solution points. 
In addition to the Radau polynomial, two other choices for g that result in new methods will be discussed. 

Next, the reconstruction for u is continuous across the interfaces, but its derivative can be 
discontinuous. For the evaluation of second derivative, we can apply the above correction procedure to the 
derivative of the reconstruction. Before this application, we need a common derivative at each interface, for 
which the simplest choice is the average derivative employed in the BR2 method. This common derivative 
was originally cast via the integral formulation; here, it is cast in a more intuitive manner via the 
differential formulation. The concept of reconstruction also clarifies and simplifies its motivation and 
calculation. Moreover, the current approach leads to several new criteria for the definition of common 
quantities. In particular, a scheme that shares the high accuracy property of Van Leer’s recovery method 
without the disadvantage of interpolating across cells is derived. The comparison and trade-offs among the 
various choices of new and existing common values as well as choices of correction functions are carried 
out by Fourier analysis. Extension to the case of quadrilateral meshes is obtained via tensor products, and 
Fourier analyses in two dimensions are discussed. For the case of the two-point boundary value problem 
(steady state), it is shown that these reconstruction schemes, which include BR2, LDG, CDG, and recovery 
(RDG), yield exact common interface quantities as well as exact cell average solutions for almost all K . 
Finally, as a side benefit of the current approach, the (continuous) reconstruction polynomials are shown, 
via numerical examples, to be more accurate than the (discontinuous) solution polynomials. 

In the hope of attracting the interest of researchers who are not familiar with these types of high-order 
schemes, this paper is written in a self-contained manner. It is organized as follows. The diffusion equation, 
solution polynomials, and various notations are introduced in Section 2. The evaluation of the first 
derivative via reconstruction is presented in §3, and that for the second derivative, in §4. Fourier analysis 
and the proof that the eigenvalues are independent of the solution points chosen can be found in §5. 
Stability and accuracy of various schemes are presented in §6. Extensions to the case of two spatial 
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dimensions are discussed in §7. Section 8 deals with numerical examples. Conclusions and discussions can 
be found in §9. Finally, the appendix concerns the boundary value problem (steady state). 


2. Diffusion Equation 

On the domain [a,b ] , consider the diffusion equation 

U t = 11 XX 


with initial condition 


z/(x,0) = u 0 (x) . 


( 2 . 1 ) 

( 2 . 2 ) 


Let the domain [a, 6] be divided into possibly nonuniform cells or elements Ej of cell widths h j and cell 
centers Xj, j . For the moment, the solution is assumed to be periodic; therefore, boundary 

conditions are trivial: the data in the ghost cells j = 0 and j =J + 1 are respectively identical to those in 
cells J and 1 . 

On each cell Ej , let the solution be approximated by the K pieces of data itj k , k =1, ..., K at locations 
Xj k , which are called solution points. The K solution points are typically the Gauss or Lobatto points. For 

(at least) linear problems, equidistant points can also be employed. In fact, it will be shown that the Fourier 
stability and accuracy results are independent of the type of points chosen. For simplicity, we assume that 
the same type and same number of points are employed for all cells. 

At each solution point, the solution u jk (t) depends on t ; Uj k (t n ), however, is abbreviated to u jk . At 


time level n, suppose the data u jk are known for all j and k. We wish to calculate diij k (t)/dt at time 
1 = t n , which is abbreviated to d u . k l dt . That is, we wish to evaluate u xx at the solution points Xj k in 

terms of the data. Then, we march in time by, say, a Runge-Kutta method. 

Instead of dealing with the cell E , or the global description, it is often more convenient to deal with the 

interval I = [-1,1] or the local description; see, e.g., (Hughes 2000). With £ varying on / and x on Ej, 
the linear function mapping I onto Ej and its inverse are 

x(g)=Xj+ghj / 2 and g(x)=2 (x-Xj)lhj. (2.3a, b) 

Let the solution points on 1 be denoted by % k , k = 1 . The solution points on Ej are, by (2.3a), 

Xj, k = *(& ) = xj + 4 hj 12. (2.4) 

A function ?y(x) on Ej results in a function on I denoted by, for simplicity of notation, r j ( c) : 

r,(£) = 0 (*(£)) = r j(x) ■ 

The derivatives in the global and local descriptions are related by the chain rule, 

dr j (*) = 2 drjtf) 
dx hj dij 

On Ej , for each k. let (/> j k be the Lagrange polynomial of degree K - 1 that takes on the value 1 at 

x jk and 0 at all other solution points. In the global and local descriptions, respectively, 

K x -x K P-P 

<t>j,k(x)= ]~J and <M£)= Y[ p _p ■ (2.6a, b) 

1=1, l*k X J’ k X J’ 1 1=1, l*k £ k 

These (j) j k are called basis (or cardinal or shape) functions. The solution polynomial, i.e., the polynomial 
of degree K- 1 interpolating u jk , k =1,..., K , can be written as, 

K K 

iij(x) = Yj u .hk ^j.kix) and u^) = ^ u j,k MS) ■ (2.7a, b) 

k= 1 k=l 

At each interface x /+ i/ 2 , let uJ +l/2 and uJ + i /2 denote the values just to its left and right, respectively, 


u j+ 1/2 = Uj(x J+ m) and «/+i/2 = u j+i(x j+ m) ■ 


(2.8a,b) 
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In general, w , +1/2 * uJ +U2 . Using the local coordinate, 

“7+1/2 = “/(D and “7+1/2 = “z+it- 1 ) • 


(2.9a, b) 


3. Evaluation of the First Derivative 

The first task is to estimate u x at the solution points Xj k . To this end, the polynomials «■ on E j 
collectively (as j varies) form a function denoted by {z/y}, which is generally discontinuous across the 
interfaces. If we ignore these discontinuities and estimate n x by (u ) x , we obtain erroneous results. Such a 
derivative estimate involves no interaction of the data between adjacent cells, and the solution never feels 
the effect of the boundary conditions. Therefore, to estimate u x , we first reconstruct n by a piecewise 

polynomial function {u^} , which is continuous at the cell interfaces and, on Ej , approximates i ij in some 
sense (the superscript ‘C’ stands for ‘continuous’ or ‘corrected’). To maintain accuracy, the polynomial z/y 

is required to be of degree K so that {u C j) x is of the same degree as z/ ; . At the solution points, (u ( j) x 
provides a derivative approximation that takes into account the data interaction. 

Common interface values. In order for {zz^} to be continuous at the interfaces, the reconstruction 

polynomials u C j and u C j +x must take on the same value at Xy +1 / 2 • Such a value is often called ‘numerical 

flux’ (see the review paper by Cockburn, Karniadakis, and Shu (2000) or Bassi and Rebay (2000)); here, 
since the property of commonality relative to the two adjacent cells is essential, and we also need a similar 
quantity for the derivative, we use the term ‘common’. Thus, at each interface, we need to define a 
common value and a common derivative. For an advection problem, between the left and right values 

“7+1/2 and “7+1/2 ’ d,e common value is typically the upwind (flux) value. For diffusion problems, we use 

the following weighted average: with 0 < k < 1 , 

“7+1/2 = KU j+m +( 1 _K ')“y+l /2 • (3-1) 

If k = 1/2 , we obtain the average (centered formula) employed by Bassi and Rebay (1997a, 2000) 

“7+1/2 = (l/2)(Wy +1 / 2 +“y+i/ 2 ) • (3.2) 

If k = 1 (or k = 0 ), we have a one-sided formula employed by Cockburn and Shu (1998), 

“/+ 1/2 _ “ 7 + 1/2 • (3-3) 

Next, we require u C j(x) to take on the values w c “” 2 at Xj_ 1/2 and u ™ ” 2 at x j+ 1/2 > t0 t> e of degree 


K , and to approximate m ; (x) in some sense (see Fig. 3.1(a)). Instead of defining z/y (x) , we define 
My (x) - zzy(x) , which approximates the zero function. Switching to the local coordinate £ on I = [-1,1] , at 
the left and right interfaces respectively, “ 7 ( 7 ) - “ 7 ( 7 ) is required to take on the following known values: 

“ 7 ( _1 ) _ “/( _1 ) = ll j°T /2 - “/( -1 ) and “7 C 1 ) — “/O) = “ 7+172 — “yC 1 ) ■ (3.4a, b) 

We now separate the prescription of the left interface value from that of the right. This separation is 
essential for the new approach to work. On the interval I = [-1,1] , let g LB be the correction function for 
the left boundary defined by the following conditions: 

(3.5) 


g LB (-l) = l, g LB (l) = 0, 


and g LB is of degree K and approximates the zero function in some sense. The condition g LB (-l) = l 


deals with the jump (3.4a), i.e., u c °™ l2 -Wy(-l) , at the left interface, while the condition g LB (l) = 0 leaves 
the right interface value z/y(l) unchanged. Let g RB be the correction function for the right boundary 
defined as the reflection of g LB : 
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«“(!) = 1, g RB (-l) = 0 > (3.6) 

and g RB is of degree K and approximates the zero function. Employing g LB , the polynomial 

«“(£) = utf) + [« 7 C T/2 -«;(-!)]*“(#) (3.7) 

provides a correction for ufE,) at the left interface by changing its value from «•(— 1) to m“” 2 , while 
leaving the right interface value w ,(1) unchanged. For later use, the polynomial 

«f(£) = + [»; + T/2 -«y(l)]g RB (l) (3.8) 

provides a correction for My(£) at the right interface by changing its value from Uj( 1) to w c °” 2 , while 
leaving the left interface value Wy(-l) unchanged. See Fig. 3.1(b). Next, the polynomial 

u C j{$) = Uj(Z) + [uJTn-Uji- l)]g LB (« + -»y(l)]g RB (^) (3.9) 

provides the corrections to both interfaces: using (3.5) and (3.6), one can easily verify that i ij (-1) = Wy“” 2 
and u C j{ 1) = 2 . See Fig. 3.1(a). Thus, the above My (£) is of degree £ , takes on the common values at 

the two interfaces, and approximates Uy(£) in the same sense that g LU and g RB approximate the zero 
function. Finally, the derivative of uj is 

(«,-)*(£) = + [z/;_T /2 -»y(-l)](g LB )'(^) + K + T /2 -«/l)](g RB )'(^) . (3.10a) 

At the solution point % k , the corrected derivative is given by 

(u x ) jk =(2lh j )(u c j)^ k ). (3.10b) 

Once the correction functions g LB and g RB are chosen on [—1,1] , the above calculations can be carried 
out in an economical and straightforward manner. 

Note that the reconstruction polynomial u'j clarifies the ideas; in practice, however, we only need 
the values of its derivative at the solution points. 



0 12 3 


(a) 



-0.5 0 0.5 1 1.5 2 2.5 


(b) 


Figure 3.1. (a) Centered common values (dots) and reconstruction functions uj and uF j for the case 
of piecewise linear solution polynomials; here, {Uj} is continuous at all interfaces x,- +1/2 - (b) Functions 
M y LB (corrected for the left boundary’) and m RB (corrected for the right boundary’). 


Correction functions. In the rest of this section, we define three correction functions. To this end, we 
need some reviews and preparations. Let the L 2 inner product of any two polynomials v and w on E ■ be 

(v,w)= ( Xj+,/2 v(f)w(f)df. 

Jx J-ll2 
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For any integer m > 0 , let P"' be the space of polynomials of degree m or less. Then P'" is a vector 
space of dimension m + 1 . When the domain, say, Ej , of the polynomials is important, we use the notation 

P"\E j ) . On / = [-1,1] , a polynomial v is orthogonal to P m = P”\I) if, for each k , 0 <k <m , 

{v,S k )=\[v{S)S k dS = o. 

Clearly, the criterion of being orthogonal to P"' provides m + 1 conditions (or equations). 

Again, on 7 = [-1,1], for £ = 0,1,2, ..., let the Legendre polynomial P k be the unique polynomial of 

degree k that is orthogonal to P k 1 and satisfies P k ( 1 ) = 1 . Thus, if k > m , then P k is orthogonal to P m . 
The Legendre polynomials are given by the following recurrence formula (see, e.g., Hildebrand 1987): 

4 = 1, 4=<T, 

and, for k > 2 , 

?£-1 k - 1 

4(£) = -V- ^4-l(£) - V Pk- 2(«- (3-11) 

k k 

The first few Legendre polynomials are plotted in Fig. 3.1(a). Properties of the Legendre polynomials that 
will be employed are listed below. First, P k is an even function (involving only even powers of L ) for 
even k , and an odd function for odd k . For all k , the values at the boundaries (or end points) are 

7>*(-l) = (-l)*, 4(1) = 1. (3.12a,b) 

The derivative values at the end points are 

7>/(-l) = (-l)*-'£(£ + l)/2, P/(l)=£(£ + l)/2. (3.13a,b) 

In addition, 

(^a ) =2/(2^ + 1) , and, for k±l, (P k ,P,) = 0. (3.14a,b) 

Finally, the zeros of P k are the k Gauss points on [-1,1]. 
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(a) (b) 

Figure 3.1. (a) Legendre polynomials and (b) right Radau polynomials. 


The right Radau polynomial of degree k ( k > 1 ) is defined by 

^=■^( 4 - 4 - 1 )- ( 3 - 15 ) 

The first few Radau polynomials are plotted in Fig. 3.1(b). Here, the letter R stands for ‘Radau’ and the 
subscript R for ‘right’; the factor (— 1)* is nonstandard and is needed so that (3.16a) below holds. 
Equation (3.15) implies that R r k is orthogonal to R k . In addition, 

and ^,-t( 1 )=°- (3.16a,b) 

Note that R r k , which is of degree k , is determined by the above two conditions and the k - 1 conditions 

that it is orthogonal to P k . This definition of the Radau polynomial shows that it is a natural choice for 
the correction function. At the two boundaries, using (3.13), we have 

*ju'(-!) = -W 4 and 7^/(1) = (-l)*-‘*/2 . 


(3.17a,b) 
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The zeros of the Radau polynomial R r k are the k right Radau points. 

The Lobatto polynomial of degree k ( k > 1 ) is defined by 

Lo a - = P k -P k - 2 - (3.18) 

Using (3.15), the Lobatto polynomial can be expressed in terms of Radau polynomials: 

= 2(-l)*(if^ >i -RR'k-i) ■ (3.19) 

The zeros of the Lobatto polynomial of degree k are the k Lobatto points; they include the two boundary 
points ±1 . As shown in Fig. 3.1(b), they are also the £ -coordinates of the intersections of the graphs of 

R R,k and R R,k- 1- 

np RR LB 

Returning to the correction functions, we always choose g by reflection: g (£) = g (-£) . As a 

result, we need to focus only on g LB . For simplicity of notation, set 

lb 

g = g ■ 

Since g is of degree K , it is determined by K + 1 conditions. Two conditions, namely, g(-l) = l and 
g(l) = 0 , are known; therefore, K - 1 additional conditions remain. We discuss below three choices for g . 

If the K - 1 additional conditions are the requirement that g is orthogonal to P K 2 , then g is the 
right Radau polynomial of degree K defined in (3.15). The resulting method yields a solution identical to 
that of DG (Huynh 2007). Therefore, we use the notation g DG : 

gDG~ R R,K- (3.20) 

The reason the DG scheme can be cast in the above reconstruction form is sketched below. In the 
current formulation, since the projections of g and g onto P (/) are zero, the projections of the 


polynomials u C j{x) (degree K ) and Uj(x) (degree K — 1) onto P k 2 (Ej) are identical. In the standard 
formulation, if (j) is a test function of degree A" - I . then <j) x is of degree K — 2 , and the expression 
| Uj(j) x dx contains the projection of iij onto P K 2 (Ej). 

Loosely put, the current formulation is a finite-difference DG formulation (versus finite-element). It 
involves no quadratures and has the advantage of simplicity and economy. In addition, regardless of the 
choice of solution points, no matrix inversion is needed. For example, if we choose the Lobatto points as 
solution points, then, since the corresponding basis functions are not orthogonal, the standard DG 
formulation requires a matrix inversion, whereas the current formulation does not. In other words, the mass 
matrix inversion is built in. 

For the next two choices, in addition to g(-l) = l and g(l) = 0, we require g to be orthogonal to 
P K ’ (yielding K-2 conditions) together with one additional condition. It was observed in (Huynh 2007) 
using Fourier analysis that the requirement of g being orthogonal to P K gives rise to stable schemes for 

convection. (The converse is not true, however: there are correction functions not orthogonal to any P k 
that result in stable schemes). 

The second choice for g requires that it vanishes at the K — \ Gauss points. (These points are the zeros 
of P K , and are completely different from the K Gauss points, which typically are the solution points). 
For obvious reason, this choice is denoted by g Ga . It will be shown later that 


gGa 


K 


2K - 1 


-R 


R,K 


K-l 
2K - 1 


R 


R,K - 1 • 


(3.21) 


Note that this g results in a scheme very similar to the staggered-grid scheme of Kopriva and Kolias 


(1996) as well as the spectral-difference method of Liu et al. (2006). The key difference is that the current 
scheme employs only one grid, and it is stable, whereas the staggered-grid and spectral-difference methods 
employ two grids, and they are mildly unstable (Huynh 2007). 

The third choice for g , requires that g' vanishes at K - I of the K Lobatto points (Legendre-Lobatto, 


not Chebyshev-Lobatto). It can be written as (compared to g Ga , the weights are switched), 
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^Lump, Lo n jr 1 Rr,K n jr . Pr.K- 1 • (3.22) 

ZK — 1 ZA — 1 

For this correction function, if the solution points % k are chosen to be the Lobatto points , then the 
derivative (i/| B )e at all solution points are the same as the uncorrected value (« ; ) f except at the left 
boundary. In other words, this g lumps the correction to the left boundary — thus, the notation g Iump; Lo . 
We will also need the following derivative values at the left boundary: by (3.17a), (3.21), and (3.22), 
Sdg'H) = -^ 2 /2, g G a'(-l) = -[^-l) + l]/2, and g Lump , Lo'(-l) = -K(K - 1)/2 (3.23) 

The plots of the three correction functions for the cases of K=2 and K = 4 are shown in Fig. 3.2. 
Loosely put, among the three correction functions, g DG is the steepest, and g Lump . Lo , least steep. Note that 
for each K , the 2 -coordinate of the intersections of the three curves are the K Lobatto points, and the 
derivative of g Lump . Lo vanishes at these points except at the left boundary <g = -\ . Also note that as K 
gets larger, g Ga gets closer to g Lump . Lo . 
















<5 Lump, Lo 







gDG^ 

'"-v^Gauss 







(a) 

Figure 3.2. Correction functions for (a) K = 2 and (b) 



(b) 

K = 4. 


Concerning practical calculations, one can derive the equations for these correction functions (for 
arbitrary K ) by a software package such as Mathematica or Maple. 

Readers who are not interested in the proofs can skip the rest of this section with no loss in continuity. 
We now discuss the derivation of the last two correction functions. Since both R rk and R R K , are 

orthogonal to P K ~ 3 , these two correction functions can be written as, with 0 < a < 1 , 

g= aR R K +l\-a)R R K _ , . (3.24) 

Due to (3.16), the above obviously satisfies the conditions g(-l) = 1 and g(l) = 0 . 

Concerning g Ga defined by (3.21), we need to prove that it vanishes at the K - 1 Gauss points. To this 
end, by expanding (3.21) in terms of the Legendre polynomials using (3.15), 

gGa - (-1)* [K Pk - (2 K - 1) P K A + (K - 1) P K _ 2 ] / [ 2 (2K - 1) ] . 

By (3 . 1 1 ), or the recurrence formula for P K , we have (2 K - 1) c P K , = K P K + (K - 1) P K 2 . Thus, 

(2K - 1) (<f - 1) P K _ X = KP k - {2K - 1) P K _ X + (K - 1) P K _ 2 . 

The above two equations imply that g Ga =[ (-1)^/2] ( c - I ) P K \ ; therefore, g Ga and ( ^ - 1 ) P K , have 
the same zeros. This completes the proof for g Ga . 

As for grump, Lo 5 b > s defined by (3.24) where a is chosen so that the point £ = 1 is a zero of 
multiplicity two, i.e., g'(\ ) = 0 . Using (3.17b), one can verify that g'( 1 ) = 0 if a = (K - 1)/(2 K - 1) . This 
a results in (3.22). It turns out that (gLump, Lo)’ vanishes not only at £ = 1, but also at K — 1 of the K 
Lobatto points. (The proof of this fact involves some algebra and can be found in (Fluynh 2007).) 
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4. Evaluation of the Second Derivative 

Unless otherwise stated, all discussions concerning stability and accuracy are based on Fourier (Von 
Neumann) analysis. Recall that for the diffusion equation, if an explicit time-stepping method is employed 
(say, Runge-Kutta), the time step size has a limit, above which the solution becomes unstable. This stability 
limit is inversely proportional to the largest magnitude of the eigenvalues of the discretization. It turns out 
that the eigenvalues are real and negative. Therefore, to compare the stability limits of schemes, we only 
need to compare the minimum eigenvalues: the scheme with a larger minimum eigenvalue (i.e., the 
minimum eigenvalue has a smaller magnitude) has the advantage of possessing a larger stability limit. 

Next, a scheme using K solution points is super-convergent (or super-accurate) if its order of accuracy 
is higher than the expected order of K . All schemes discussed here are super-convergent for K> 3 : their 
orders of accuracy are at least 2 (K - 1) . In general, super-convergence does not carry over to nonlinear 
equations: due to nonlinear errors, the order of accuracy is typically only K . For comparison, we discuss 
the main trade-offs in accuracy and stability here; additional details will be provided in §6. In general, a 
less steep correction function for the solution points results in a scheme with the advantage of a larger 
stability limit and the disadvantage of lower accuracy. 

Returning to our current task, recall that {i/j \ is continuous across the interfaces; as a result, {u C j) x 
provides an approximation to u x that takes into account the data interaction. For each cell Ej, at the 
solution points x jk , set 

v j,k=( u j)x( x j,k)- (4-1) 

These K (derivative) values define a polynomial of degree K — 1 denoted by Vj{x) that is identical to 

(u < j ) x (x) . Note that the function {v } = {v y (x)} = {( it j) x } can be and usually is discontinuous across the 
interfaces. Therefore, to calculate u xx , we can apply the reconstruction procedure once again to the 
discontinuous function {v ; J = {{Uj) x } . To this end, we must define a common derivative (or common 
gradient) at the cell interfaces. 

Common derivatives. At each interface, in formula (3.1) for the common value, with 0 <k<\, the 
weight for u ~ +1/2 is k and that for u + +l/2 , \ — k. To define the common derivative, we switch the two 

weights. Loosely put, this switch serves the purpose of canceling out the one-sidedness so that the scheme 
is consistent with the centered or unbiased nature of the diffusion process. Since the corrected derivative 
(; u C j ) x is readily available, an obvious choice is 

( U x)j+ 1/2 = 0 ~ v)(m / ) x (Xj + i/ 2 ) + K ( u j+l)x( x j+V2 ) • (4.2) 

The function Uj involves the data in the cells j — l, j , and / + 1 . Consequently, the above formula has a 
four-cell stencil involving the data from E: , to Ej +2 as can be seen in Fig. 4.1(a) for the centered case 
( k = 1/2). Since the calculation of in the cell Ej employs (u x )j°y 2 and (t/ v )““ 2 , the corresponding 
scheme has a five-cell stencil. 

We now define a common derivative at each interface x J+1/2 in a manner that it involves only the data 
in the two adjacent cells, namely, Ej and E . A scheme with such a compact stencil is desirable since it 
is easy to code, the boundary conditions involved are simple, and the resulting implicit version has a 
sparser as well as invertible matrix. With i/y®(x) by (3.7) and uJ H (x) by (3.8), set 

K )jlv2 = (l-*)(«f ),(*, +1/2 ) + *(« “)*(*, + l/2) ■ (4-3) 

The functions n^ B and uj +i for the centered case are shown in Fig. 4.1(b) and, for the one-sided case with 
k - 1 , Fig. 4.2(a). The above can be expressed as: 
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0',)“w2 = (!-«•) f t («>)#(!)+ [« y T/2 "«,(!)] [(*“)' (1)] 1 + 


7+1 


(«y +1 )^(-l)+ K + T/ 2 -^ +1 (-l)] fe LB ) (-!)]/■ 


(4.4) 



0123 0123 


(a) 


(b) 


X 


Figure 4.1. Centered common value and derivative, (a) = (1/2) [(Wy) x (lC y+ i/2) + ( u< j+l) x ( X j+\l2)\ 

has a four-cell stencil, (b) (u x )Jffi 2 = (1/2) [(m^)* 07+1/2) + (t7.f1)* (-*7+1/2)] f Jas a two-cell stencil. In 
both figures, the solution polynomials are linear, and the correction function is g DG . 


for the definition of the common interface 

is to be avoided 


As for the boundary corrections (g LB ) (-1) and (g RB ) (1) , they can be expressed in a formula. Recall 
that due to symmetry, (g RB ) (l) = -(g LB ) (-1) . The three correction functions in the order of increasing 
steepness are g Lump ^ , g Gauss , and g DG . The corresponding (g^Ml) are, respectively, by (3.23), 

K(K- 1)/2, [K(K- l) + l]/2, and K 2 / 2. 

Note that for steady state problems, if we use g Lump Lo 
derivative, the resulting matrix for the solution of u j k becomes singular. Thus, g L Lo 

for steady state case. (It can still be used for the correction at solution points, however.) 

Also note that at the common interface x J+l / 2 , the common derivative evaluated by (4.3) or (4.4) 

involve corrections only at that interface. On the other hand, at the solution points of Ej , the derivative 
evaluated by (4. 1) using u - involves corrections at both interfaces of E • . 

With the common derivative defined by (4.4), set 

com / \ com 

V j+H2 “ ( U x)j+ 1/2 ' 

Finally, by applying the reconstruction procedure to the data v,- k of (4.1) and the above common 


quantities, at each solution point, the corrected second derivative is given by 

( u xx)j,k=(. v j)x( x j,k) ■ (4.5) 

A few remarks are in order concerning the trade-offs between (4.2), which has a wide stencil, and (4.3), 
which has a compact stencil. 

If the centered common value is employed, then as shown in Fig. 4.1, expressions (4.2) and (4.3) yield 
different schemes. These two schemes are comparable in terms of accuracy and stability; the scheme using 
(4.2) is slightly more accurate but its stability limit is smaller. However, the disadvantages of (4.2) are 
numerous: its stencil is large and thus, the coding of boundary conditions is more complex; it has a zero 
spurious eigenvalue with an undamped eigenfunction; and the steady state case results in a singular matrix. 
Therefore, unless otherwise stated, we employ only the compact stencil formula, namely, (4.3). 

If g DG is employed as correction functions at both the solution points and interfaces, then, with the 
centered common values, (4.2) results in a method identical to BR1, and (4.3), BR2 (Bassi and Rebay 
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1997a, 2000). These two BR schemes were originally formulated in the finite -element framework whereas 
the current versions are formulated in the finite-difference one. 

In the case of common quantities via one-sided choice (e.g., the value to the left and the derivative to 
the right), for the one -dimensional case discussed here, the schemes corresponding to (4.2) and (4.3) are 
identical. For the two-dimensional case, however, they can be different. If, in addition, g DG is employed as 
correction functions for both the interfaces and the solution points, then (4.2) results in a method identical 
to the local discontinuous Galerkin or LDG (Cockburn and Shu 1998), and (4.3), the compact 
discontinuous Galerkin or CDG (Peraire and Persson 2008). Again, these DG schemes were originally 
formulated in the finite-element framework whereas the current versions are formulated in the finite- 
difference one. 

Comparing the centered versus one-sided common values with g DG as correction function (i.e., BR2 
versus CDG), the former is of order 2 (K -1) ; the latter, order 2 K ; the former, however, has the advantage 
that its stability limit is more than two times larger than the latter. Since super-convergence (accuracy 
higher than K ) does not hold for the general case of nonlinear equations, the scheme using centered 
common values appear to have an edge. (For additional comparisons, see the numerical examples in §8.) 

Criteria with high accuracy for common derivatives. Continuing with the derivation of the common 
derivative, new criteria with the advantage of high accuracy can be obtained as follows. 

We reconstruct the solution on the domain Y x j-\n^ x j^sn\ — Ej U£ y+ i in a manner that (a) it takes on 
the common value which is the weighted average (3.1) at x- +1/2 , and (b) it approximates the data on E j 
and E +l without the requirement of leaving the interface values at x f l/2 and x +3/2 , respectively, 
unchanged. In other words, for the cell E +] , to correct for its left boundary x - +1/2 , with g = g LB , we 
require that g(-l) = 1 , while g(l) is allowed to be arbitrary and, critically, we require g to approximate 
the zero function to a degree as high as possible. A similar statement holds for the cell E ■ and its right 
boundary x y - +1/2 . Thus, to calculate the common derivative (but not the derivative at the solution points) 

we can employ the Legendre polynomial of degree K , which is orthogonal to P K , as the correction 
function: 

gu=g LB =(-l)*Pjr and g RB =^. (4.6) 

By (3.13b), for this connection function, 

(^ RB )'(1) = -0f LB )'(-l) = K(K + 1)/2. 

The common derivative is again given by (4.4). 

For the case of centered common values using g DG as correction function at the solution points, the 
method using (4.6) has the following order of accuracy: for odd K , the order is 2 K , and for even K , 
2 (K - 1) ; e.g., both the schemes with K = 3 and K = 4 are of order 6. 



0 12 3 



0 12 3 


(a) (b) 

Figure 4.2. (a) The common derivative (w x )y°” 2 = (z/ BB j) x (x y+1/2 ) , which corresponds to (4.3) with k = 1, 

has a two-cell stencil, (b) Using Van Leer and Nomura ’s recovery, the common value ii(xy +1/2 ) an d 
common derivative w x (x y+1 / 2 ) have a tw’o-cell stencil. 
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We present now a new criterion for the common quantities (both value and derivative) that is centered 
(or unbiased) and, with g DG as correction function at the solution points, yields a scheme of order 2 K for 

all K . At each interface Xj +U2 , let z/““ 2 be an unknown. Using the Legendre polynomial (4.6) as 
correction function for the interface (a different g results in a lower order of accuracy), the corresponding 
corrected derivatives (w^ B ) JC (jc . + 1/2 ) and {Uj +l ) x {Xj +ln ) can easily be obtained as linear functions of 
Uj+T/ 2 ■ The unknown « 2 is then determined by requiring that these two derivatives are equal: 

( 11 j )x( x j+iii)- (. u j+i)x( x j+m) ■ (4-7) 

In other words, on (Xj_ l/2 ,Xj +3/2 ) > both the reconstruction and its derivative are required to be continuous 
across x - +1/2 . (Note that the continuity requirement simultaneously for all interfaces involving a 

tridiagonal solver was briefly discussed by Kopriva (1998) for his staggered-mesh methods. The current 
method is considerably simpler, however.) 

For our final criterion, we define both Wy““ 2 and (u x )J° t y 2 via the recovery approach of Van Leer and 

Nomura (2006). Let the recovery function u be a polynomial of degree 2 K - 1 on the domain 
[ x j-U2’ x j+3n] = Ej U^y+i determined by the 2 K conditions that: it and w ■ have the same projection on 

P k ~ l (E ■) , and it and a +1 have the same projection on P lc ~ 1 (Ej +1 ) . In other words, for 0 < i < K - 1 , 

J ii(x) (x — x ■)' dx = u Ax)(x —x ■)' dx, (4.8a) 

E : j 

and 

f u(x) (x - x +1 ) ! dx = [ u i+l (x)(x-x j+1 )' dx . (4.8b) 

The common quantities at the interface are then defined as u(x j ];2 ) and (u x )(x , + i/ 2 ) ■ See Fig. 4.2(b). 

Compared to the averages (3.1) and (4.3), the common quantities via recovery above has the advantage 
of higher accuracy and the disadvantage of being costly. In fact, based on Fourier analysis for K up to 10 
(with K > 10 , the errors become exceedingly small), the gain in accuracy is considerable: using g DG as 
correction function at the solution points, for odd K , the order of the recovery scheme is 3 K - 1 , and for 
even K , 3K-2 . A peculiar property of these methods is the following: for K> 4 , two of the spurious 
eigenvalues become complex for a certain range of wave numbers, which means that the corresponding 
components are propagated as well as damped. These complex eigenvalues, however, pose no problem 
since the corresponding components are damped very quickly. Also note that we can consider the averages 
(3.1) and (4.3) as lower order approximations to u(x /+l/2 ) and (it x )(x J+1/2 ) , respectively. 

Algorithm. The algorithm for the diffusion equation is summarized below. Let the data itj k at the 
solution points Xj k be given for all cells. 

fit At each interface Xj +U2 , calculate uJ +1/2 and uJ +l/2 by (2.8a, b) and the common value 

(3.1). Next, obtain the common derivative (u x )j7i/2 by (4.4). Alternatively, «“™ 2 and {u x )f™ 2 
can be estimated via the continuity requirement (4.7) or via the recovery it defined by (4.8). 

□ For each cell, at the solution points, calculate Vj(x jk ) = (u C j) x (Xj k ) via (3.10). With 

v 5 +T /2 = ( u x)j+ T /2 at the interfaces, the values of u xx at the solution points are given by 

(v C j)x( x j,k) ■ 


Boundary conditions. We end this section with a discussion of boundary conditions. Let the left 
boundary of the domain be x /+1/2 for some / . If we have a Dirichlet condition u D at this boundary, then 

we define the common value by 


13 


com - / a r\ \ 

u i+ 1/2 _ u 1+112 ~ u d ■ (4.9a) 

This definition corresponds to k = 1 in (3.1). Next, using (4.3) for the common derivative, set 

(»XT 2 = (0,(*/ +1 /2). (4.9b) 

On the other hand, if we have a Neumann boundary condition u x = ( u x ) N at the left boundary jc /+ 1/2 for 

some real number (u x ) N , then we define 

( ll x ) l + U 2 ~( ll x ) l + 1/2 ~( ii x)n • (4.10a) 

The value M/ C °y 2 * s determined by the condition that («/+i ) x i x i+\n) ~ % > *- e -> 

u i+ 1/2 = w /+i(~l) + P, + i/2)« w -(M /+1 ) # (-l)]/(g ) ( 1) • (4.10b) 

The boundary conditions on the right end of the domain follow by reflection. 

As will be shown later, for the case where the exact solution is a parabola, and two Gauss points (linear 
element) are employed in each cell, the standard DG schemes do not recover the exact solution. Together 
with the above boundary conditions, the two schemes that recover the exact solution are: the scheme using 
the Legendre polynomial with the continuity requirement (4.7) and the recovery scheme (4.8). 


5. Fourier (Von Neumann) Analysis 

Fourier analysis provides information on both stability and accuracy. On (-oo, oo), consider the equation 
(2.1), i.e., Uf — u xx , with the periodic initial condition w init (x) = e lwx , where the wave number w lies 
between -n and n . Low frequency data corresponds to w of small magnitude; high frequency, to w 
near ±n . The exact solution is M exact (x,l) = e~" ‘e ,wx . At (x,t) = (0,0) , 

(»exac,),(0,0) = -H’ 2 . (5-1) 

The cells are Ej = [j - 1 12, j + 1 / 2] . The data are it j k = exp [ i w(j + / 2) ] . However, it is not the data 

but their following property, which plays a key role in the calculation of eigenvalues, 

11 j-\,k ~ e lw i l j, k ■ (5.2) 

To calculate the eigenvalues, the K solution points are grouped together as a vector: set 


“j = 


(5.3) 


li : is 

v 


For the case of three-cell stencil, the solution can be expressed as 

d ll j ^ ,. 


(5.4) 


Using (5.2), we replace Uj +l by e llw u^ , 


du 

dt 


L = { J j e^C l )u j 


i=-\ 


where C) are K ■ K matrixes. Thus, the spatial discretization (semidiscretization) results in 

— — = S iij where S = V e'^’Cj . 
dt 1 1 


(5.5) 


(5.6, 5.7) 


Here, S, for ‘space’ or ‘semidiscrete’, is a K ■ K matrix with K eigenvalues. The one approximating - vv 2 
(the exact eigenvalue for d 2 / dx 2 ) is called the principal eigenvalue and is denoted by S(w) : 

S(w)«-w 2 . (5.8) 

All others are spurious eigenvalues. For stability, all eigenvalues must have nonpositive real parts. 

If an explicit time-stepping method such as the standard four-stage Runge-Kutta is employed, then the 
stability limit of a scheme is proportional to the inverse of the largest magnitude of all eigenvalues. 
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Accuracy. Concerning accuracy, note first that for a uniform mesh of cell width h , S approximates 
h 2 d 2 Idx 2 . A scheme is accurate to order m if S approximates h 2 d 2 /dx 2 to 0(h m+2 ) , i.e., for small w, 

S(w) = —w 2 + 0(w m+2 ) . (5-9) 

In practice, it is difficult if not impossible to derive an expression for the principal eigenvalue A(w) when 
the number of solution points is greater than 2 or 3. Therefore, we obtain the order of accuracy by the 
following approximation. Let w be fixed, say, w = /z78 . We calculate the error 

E{yv) = S(w) + w 2 . (5.10) 

By halving w (i.e., doubling the number of cells), the error corresponding to w/2 (say, w/2 = 7r/16)is 

£(w/2) = A(w/2) + (w/2) 2 . 

Since 0((w/ 2)" i+2 ) = (1/2 ) m+2 0(w m+2 ) , for a scheme to be / 77 -th order accurate, 

E(w/2) *(l/2)" ,+2 £(w) . (5.11) 


That is, E(w)/E(wl 2) ~ 2 m+ . Using the logarithm function, the order of accuracy is given by the integer 


Log 

Lf"> Uogffl 


{E(w/2)J/ 


-2 . 


(5.12) 


Choice of solution points. Next, we discuss the choice of solution points. The two common choices are 
the Gauss and Lobatto points. While these choices may be advantageous for the nonlinear case, from the 
consideration of linear stability and accuracy, all choices are equivalent. In fact, we claim that the 
eigenvalues of A are independent of the solution points chosen. 

To prove this claim, consider the local coordinate E, on I =[-1,1] with solution points E k and data u k , 

k = l,...,K . Denote the vector {u k } k=x by u. Recall that the solution polynomial is p = yf u k <j> k {^) . 

Next, let , l = 1, ..., K be another set of solution points (we assume no two points are the same), and let 

Uj be the value p(Ei) , 

Set u = \uj}f = j . In addition, for 1 <k,l < K , set 

»*/* =&(?/)• (5.14) 

Denote the K ■ K matrix \m lk \ by M . Then, it follows from the above two expressions that 

u=Mu. (5.15) 

Since we can obtain u from u , M is invertible. Next, equation (5.7a) takes the form, 

du/dt=Su. (5.16) 

For the solution points D , the corresponding equation is 

dill dt = S u . (5.17) 

We wish to show that S and S are similar, and thus, have the same eigenvalues. To this end, multiplying 
(5.16) on the left by M , we have 

M du / dt = MS 11 . 

Since M is independent of t and is invertible 

d(M u)l dt = MSM X M u = (MSM~ l )(M u) . (5.18) 

Therefore, by (5.18), 

du / dt = (MSM- X )u . (5.19) 

Comparing (5.17) with the above and, since « is arbitrary, we have 

S = MSM 1 . (5.20) 


This completes the proof. 
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6. Stability and Accuracy of Schemes for Diffusion 

In this section, we employ the Fourier analysis of the previous section; therefore, the mesh is uniform, 
and the cell width is 1. With K solution points, all schemes discussed have order of accuracy of at least 
2 (K - 1) . As a consequence, for K > 2 , the order of accuracy is at least 2. For K = 1 , however, the order 
of accuracy can be 0, which implies that the corresponding scheme is inconsistent. Again, concerning the 
super-convergence or super-accuracy property (order higher than K ), for the nonlinear case such as the 
Navier-Stokes equations, due to nonlinear errors, it is expected that the schemes discussed are accurate to 
order no higher than K . Unless otherwise stated, all schemes have real negative eigenvalues ( < 0 ). 

A scheme is determined by the following choices. First, for the common quantities at each interface, the 
four choices are: (a) centered, (b) one-sided, (c) continuous (condition (4.7) or the derivative values to the 
left and right of the interface being equal), and (d) recovery, namely (4.8). Next, for the cases (a), (b), and 
(c), the derivative at the interface requires a correction function; in the order of decreasing steepness, the 
choices are: g he , g DG = g Ra , g Ga , and g LumPi Lo- To denote the scheme using a certain choice at the 

interface, we use the following self-explanatory notation: I - centered - g Ga , or I - continuous - g DG , or 
I - recovery . Finally, for the derivative values at the solution points, the correction functions in the order of 
decreasing steepness are: g DG , g Ga , and g Lump Lo . The corresponding schemes are denoted by, SP - g DG , 
etc... We refer to a scheme by the choice at the interface and that at the solution points, e.g., 

I- centered -g DG /SP-g Ga . 

That is, the common quantities are centered, and the derivatives just to the left and right of the interface are 
calculated via g DG ; for the solution points, the correction function is g Ga . 

6.1. The case K- 1 . Here, 

g — g I )G g(xa gLump, Lo — ^ . 

Therefore, there is only one correction function for the solution point. At the interface, however, there are 
two correction functions: the above g and the Legendre polynomial g Le = -E, . 

The three schemes I - centered - g Le , I - recovery , and I - one - sided - g DG yield a result identical to the 
standard second-order accurate method 

(.Uxx)j=U j+l +U j _ l -2ll j . 

Next, for the scheme I - centered - g DG / SP - g DG , i.e., BR2 scheme, we have u™™ /2 = (Uj +it J+1 )/2 . 
Therefore, 

K ) J + i /2 = «/+ 1 “ m 7 + T /2 = («/+ 1 - «/)/ 2 and OJ /+ 1/2 = M “?/ 2 - «/ = Oy+i “ 11 yV 2 • 

Thus, (u x ) c p , 2 = (m /+ , - Uj) 12 , and 

( u xx)j = ( u j+i +u j-i -2iij)/2 . 

This quantity approximates u xx 12 . As a result, this scheme has an error of 0(1) and is inconsistent. 

6.2. The case K = 2 . All correction functions are different from each other (see Fig. 3.2(a)). The plots 
of eigenvalues as functions of wave numbers for six representative schemes are shown in Fig. 6.1. 

The minimum eigenvalues, orders of accuracy, and errors of 18 schemes (with K = 2) are shown in 
Table 6.1. The last column shows the value g'(-l) where g is the correction function for the left interface. 
Note that all schemes are of order 2 or higher. The three schemes 11, 14, and 15 are fourth-order accurate. 
Schemes 4, 7, and 10, which employ g Lump Lo at the interfaces, have a zero spurious eigenvalue at w = 0 ; 

thus, the corresponding eigenvector is not damped. (For the steady state case, such a scheme results in a 
singular matrix for the solution and is to be avoided.) Scheme 10 is a limiting case: the spurious and 
principal eigenvalues become identical for all w . Scheme 15 (which is identical to CDG and LDG) has the 
drawback that the stability limit is rather small (the magnitude of the minimum eigenvalue is second largest 
on the list). Scheme 0 or I -centered - 3g DG /SP - g DG has the smallest stability limit of all with no 
improvement in accuracy. In general, a steeper correction for the derivative at interface results in a scheme 
with a smaller stability limit. From here on, therefore, we will not use a correction function steeper than 
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g Le at the interface. As a limiting case, scheme I -continuous- g Lump Lo /SP -g Lump -Lo (not listed) is 
inconsistent: it has an eigenvalue that is identically 0 for all w . 




(a) I - centered - g^/SP - g DG (b) I- centered -g DG /SP-g DG (c) I - centered - g Ga /SP - g Lump Lo 





(d) I - continuous - gj^/SP - g DG (e) I - recovery/SP - g DG I-one-sided-g DG /SP-g DG 

Figure 6.1. Eigenvalues as functions of wave numbers for six representative schemes with K = 2 . 


Scheme 
for K = 2 

Min. 

Eigval. 

Ord. 

Acc. 

Error 

For I-g, 
g'(-l) 

0. I- centered -3g DG /SP-g DG 

-60. 

2 

- w 4 /12 + 0 ( w 6 ) 

-6 

1. I- centered -g Le /SP-g DG 

-24. 

2 

- w 4 / 12 + 0(w 6 ) 

-3 

2. I- centered -g DG /SP-g DG 

-13. 

2 

— w 4 /12 + 0 ( w 6 ) 

-2 

3. I- centered -g Ga /SP-g DG 

- 12 . 

2 

- w’ 4 /12 + 0(w 6 ) 

-1.5 

4.1- centered - g Lump Lo /SP - g DG * 

- 12 . 

2 

w 4 /12 + 0 ( w 6 ) 

-1 

5. I- centered -g DG /SP-g Ga 

-9.7 

2 

- w 4 /24 + 0(w 6 ) 

-2 

6. I- centered -g Ga /SP-g Ga 

-8. 

2 

- u' 4 /24 + 0(w 6 ) 

-1.5 

7. I - centered - g Lump Lo /SP - g Ga * 

-8. 

2 

w 4 /12 + 0(w 6 ) 

-1 

8. I- centered -g DG /SP-g LumpjLo 

-8. 

2 

w 4 / 12 + 0(w 6 ) 

-2 

9. I- centered -g Ga /SP-g LumpLo 

-6. 

2 

w 4 /12 + 0(w 6 ) 

-1.5 

10. I - centered - g Lump Lo /SP - g Lump , ^ * 

-4. 

2 

w 4 / 12 + 0{w 6 ) 

-1 

11. I- continuous -g Le /SP-g DG 

-24. 

4 

w 6 /1440 + O(w 8 ) 

-3 

12. I - continuous - g DG /SP - g DG 

- 12 . 

2 

w 4 / 24 + 0(w 6 ) 

-2 

13.1- continuous - g Le /SP - g Ga 

-16. 

2 

w 4 / 24 + 0(w 6 ) 

-3 

14. I-recovery/SP-g DG 

-15. 

4 

w 6 / 360 + 0(w 8 ) 

* 

15. I - one - sided - g DG /SP - g DG 

-36. 

4 

w 6 /540 + O(w 8 ) 

-2 

1 6. I - one - sided - g Lump Lo /SP - g Lump Lo 

-10.5 

2 

w 4 H> + 0(w 6 ) 

-1 

17. I - centered - g DG (4 - cell)/SP - g DG 

-16. 

2 

-w 4 ! 24 + 0(w 6 ) 

-2 


Table 6.1. Minimum eigenvalues, orders of accuracy, and errors for schemes with K — 2 . 
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Figure 6.2. Eigenvalues as functions of wave numbers for six representative schemes with K = 3 . 


Scheme 

Min. 

Ord. 

Coarse mesh 

Fine 

mesh 

For I-g, 

for K = 3 

Eigval. 

Acc. 

error, 
w-n! 8 

error, 
w = 7r/l6 

g'(-l) 

1. I- centered -g Le /SP-g DG 

-75. 

6 

-5.06 

• 10~ 8 

-1.97 

• io- 10 

-6 

2. I- centered -g DG /SP-g DG 

-60. 

4 

2.15- 

10^ 6 

3.40 

10~ 8 

-4.5 

3. I- centered -g Ga /SP-g DG 

-60. 

4 

5.14- 

10^ 6 

7.97 

10~ 8 

-3.5 

4.1- centered - g Lump Lo /SP - g DG * 

-60. 

4 

7.8- 

lO^ 6 

1.2- 

10~ 7 

-3 

5. I- centered -g DG /SP-g Ga 

-42. 

4 

5.6- 

10~ 6 

8.72 

10~ 8 

-4.5 

6. I- centered -g Ga /SP-g Ga 

-36. 

4 

8.62- 

10^ 6 

1.33 

10~ 7 

-3.5 

7. I - centered - g Lump> Lo /SP - g Ga * 

-36. 

4 

1.13- 

10~ 5 

1.73 

10~ 7 

-3 

8. I- centered -g DG /SP-g LumpLo 

-42. 

4 

9.93- 

10~ 6 

1.54 

10~ 7 

-4.5 

9. I- centered -g Ga /SP-g LumpLo 

-30. 

4 

1.3- 

10~ 5 

2.- 

10~ 7 

-3.5 

10. I - centered - g Lump Lo /SP - g Lump Lo * 

-26. 

4 

1.57- 

10~ 5 

2.4- 

10~ 7 

-3 

11. I- continuous -g Le /SP-g DG 

-74. 

6 

1.87- 

10^ 9 

7.31- 

10- 12 

-6 

12. I - continuous - g DG /SP - g DG 

-42. 

4 

2.18- 

10~ 6 

3.41 

10~ 8 

-4.5 

13. I - continuous - g Le /SP - g Ga 

-60. 

4 

3.4- 

10~ 6 

5.31 

10” 8 

-6 

14. I - recovery/SP - g DG 

-33. 

8 

4.75- 

10 -11 

4.64- 

10~ 14 

* 

15. I - one - sided - g DG /SP - g DG 

-148. 

6 

4.5- 

10~ 9 

1.75- 

10 -11 

-4.5 

16. I - one - sided - g Lump , Lo /SP - g LumPi Lo 

-76.5 

4 

1.55- 

10~ 5 

2.39 

io- 7 

-3 

17. I - centered - g DG (4 - cell)/SP - g DG 

-65. 

4 

6.84 

10~ 8 

2.64- 

io- 10 

-4.5 


Table 6.2. Minimum eigenvalues, orders of accuracy, and errors for schemes with K = 3 . 


18 


0 

-25 

-50 

-75 

-100 

-125 

-150 

-175 


(a) I- centered -g Le /SP-g DG 




(b) I- centered -g DG /SP-g DG 



(c) I - centered - g Ga /SP - g Lump , Lo 


0 

-25 

-50 

-75 

-100 

-125 

-150 

-175 


(d) I- continuous -g Le /SP-g DG 




(e) I - recovery/SP - g DG 



(f) I - one - sided - g DG /SP - g DG 


Figure 6.3. Eigenvalues as functions of wave numbers for six representative schemes with K = 4 . 


Scheme 
for K = 4 

Min. 

Eigval. 

Ord. 

Acc. 

Coarse mesh 
error, 

w=nl 8 

Fine mesh 
error, 
w = n ! 16 

g'(-l) 
(I -g) 

1. I- centered -g Le /SP-g DG 

-187. 

6 

-5.31- 10 -9 

-2.16 

■ 10 11 

-10 

2. I- centered -g DG /SP-g DG 

-170. 

6 

-5.06- 10~ 9 

-2.14 

10“ u 

-8 

3. I- centered -g Ga /SP-g DG 

-170. 

6 

-3.78- 10^ 9 

-1.98 

10~ n 

-6.5 

4.1- centered - g Lump Lo /SP - g DG * 

-170. 

6 

4.5- 10~ 9 

1.75- 

10 11 

-6 

5. I- centered -g DG /SP-g Ga 

-122. 

6 

-9.93- 10~ 10 

-5.05 

10~ 12 

-8 

6. I- centered -g Ga /SP-g Ga 

-98. 

6 

below 

-3.89 

io- 12 

-6.5 

7. I - centered - g Lump Lo /SP - g Ga * 

-98. 

6 

6.64- 10~ 9 

2.59- 

10~ n 

-6 

8. I - centered -g DG /SP-g Lump Lo 

-122. 

6 

2.19- 10~ 9 

7.63- 

10~ 12 

-8 

9. I - centered -g Ga /SP-g Lump Lo 

-89. 

6 

2.98- 10~ 9 

8.55- 

10- 12 

-6.5 

10. I - centered - g Lump Lo /SP - g Lump> Lo * 

-83. 

6 

8.43- 10~ 9 

3.29- 

10~ u 

-6 

11. I- continuous -g Le /SP-g DG 

-183. 

8 

8.57- 10~ 13 

8.38- 

10" 16 

-10 

12. I - continuous - g DG /SP - g DG 

-122. 

6 

2.24- 10~ 9 

8.76- 

10~ 12 

-8 

13.1- continuous - g Le /SP - g Ga 

-170. 

6 

4.21- 10~ 9 

1.64- 

10 11 

-10 

14. I - recovery/SP - g DG 

-68. 

10 

-1.77- 10~ 14 

-4.37 

io- 18 

* 

15. I - one - sided - g DG /SP - g DG 

-439. 

8 

1.96- 10~ 12 

1.92- 

io- 15 

-8 

16. I - one - sided - g Lump> Lo /SP - g Lump> Lo 

-272. 

6 

1.50- 10~ s 

5.85- 

10^ n 

-6 

17. I - centered - g DG (4 - cell)/SP - g DG 

-176. 

6 

-1.39- 10~ 9 

-5.46 

• 10~ 12 

-8 


Table 6.3. Minimum eigenvalues, orders of accuracy, and errors for schemes with K = 4 . 

6.3. The case K = 3 . The plots of eigenvalues as functions of wave numbers for six representative 
schemes are shown in Fig. 6.2. 

The minimum eigenvalues, orders of accuracy, and errors of 17 schemes for K = 3 are shown in Table 
6.2. Note that all schemes are of order at least 4 or 2 (K -1) . The three schemes 1, 11, and 15 are of order 
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6, and scheme 14 (recovery), order 8. Schemes 4, 7, and 10, which employ g LumpLo at interfaces, have a 
zero spurious eigenvalue at w = n (not at vv = 0 as in the case of K = 2 ); the corresponding eigenvector 
is not damped; and the matrix for steady state equation is singular. Schemes I - continuous - g Lump Lo (not 
listed), regardless of the choice of SP - g , are inconsistent: they have an eigenvalue that is identically 0 for 
all w . The recovery scheme (I - recovery/SP - g DG or scheme 14) is more complicated than the others but, 
among the schemes of Table 6.2, it has the advantages of highest accuracy and third largest CFL limit (the 
amplitude of the minimum eigenvalue is 33). 

6.4. The case K = 4. The plots of eigenvalues as functions of wave numbers for six representative 
schemes are shown in Fig. 6.3. 

The minimum eigenvalues, orders of accuracy, and errors of 17 schemes for K = 4 are shown in Table 
6.3. Again, the last column shows the value g'(-T) . Note that all schemes are of order at least 6 or 
2{K-\) . Schemes 11 and 15 are of order 8, and the scheme 14 (recovery), order 10. Schemes 4, 7, and 10, 
which employ g Lump Lo at interfaces, have a zero spurious eigenvalue at w = 0 ; the corresponding 
eigenvector is not damped; and the matrix for steady state equation is singular. For scheme 6, or 
I - centered - g Ga /SP - g Ga , the error crosses zero near w-k! 8. Therefore, to calculate the order of 
accuracy, we set w = n !\ 6 for the coarse mesh and w = n / 3 2 for the fine mesh. The errors are 
respectively -3.89- 10 -12 and -1.98- 10 14 , and the scheme is of order 6. Schemes 

I - continuous -g Lump Lo (not listed), regardless of the choice of g for the solution points, are inconsistent: 

they have an eigenvalue that is identically 0 for all w . Note that for a certain range of w , two spurious 
eigenvalues of the recovery scheme become complex — a rather peculiar property. Since the real parts of 
these complex values are roughly -45 , the corresponding eigenvectors are damped quickly. The property of 
possessing complex eigenvalues is unique for the recovery scheme and it holds true for all K > 4 (this 
author tested for K up to 10 ). 

6.5. Arbitrary K . For a majority of schemes, the author tested with values K of up to 10. The 
following conclusions can be drawn. 

All schemes are of order at least 2( K - 1) . In the order of increasing accuracy, scheme 1 or 
I - centered - g Le /SP - g DG is of order 2 K for odd K and 2(K - 1) for even K . Schemes 11 and 15, i.e., 
I - continuous - g Le /SP - g DG and I - one - sided - g DG /SP - g DG , respectively, are of order 2 K for all K . 
Scheme 14 or I - recovery/SP - g DG is of order 3 K - 1 for odd K and 3 K - 2 for even K . 

In general, g Lump Lo is a limiting case in the sense that (a) a less steep correction function may lead to 
instability, and (b) for the steady state case, the scheme I - g Lump Lo results in a singular solution matrix. 

Based on the above results using Fourier analyses together with the fact that the super-convergence 
(super-accuracy) property does not hold for the general case of nonlinear systems, the following 
recommendations can be made. Scheme 2 (or I - centered - g DG /SP - g DG or BR2) and 15 (or 
I - one - sided - g DG /SP - g DG or CDG) appear to be good choices: they are simple to code and, concerning 
the stability limit, the BR2 method has an edge. Scheme 9 or I -centered - g Ga /SP - g Lllmp jLo is simple to 
code and, while not quite as accurate as the above two, has a rather large stability limit. If accuracy is a key 
concern, then at the cost of some additional coding, scheme 11 or I - continuous - g Le /SP - g DG is of order 
2 K . Among all methods, scheme 14 or I - recovery/SP - g DG is the most accurate and has the largest 
stability limit for K > 4 , but it is also the most complicated since the recovery is across two cells (not one 
at a time as the other methods). 


7. Two Dimensional Extensions 

For the two-dimensional (2D) extension to a quadrilateral mesh of the current reconstruction approach, 
the solution procedure reduces to a series of one-dimensional (ID) operations in a manner similar to the 
extension of the staggered-grid scheme in (Kopriva 1998). The key difference is that the current extension 
is simpler since it involves only one grid instead of three. 

With u - u{x, y, t ) , consider the 2D diffusion equation 
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“/=“«+«; KV (7.1) 

with initial condition 

z/(x,j,0) = z/ Mt (x,y) . (7.2) 

For the case of the Navier-Stokes equations, we have 

u t + V. F( u, Vw) = 0 

where « is the (column) vector of conservative variables and V = F( w,Vw) is the flux vector. Next, 
F = F u { a) + F v ( u. V// ) where and F v are the advective and viscous flux vectors respectively. In 
addition, F =(f,g) where / and g are the two Cartesian components; similarly, F a =(f a ,g a ) and 
F v = ( f v ,g v ) • The above equation can be written as 

U, + /v + 8y = 0 • (7.3) 

Clearly, (7.1) is a special case of (7.3): with u = m and F( u,Vu ) = Vm , i.e., / = u x and g = t/ v , equation 
(7.3) takes the form 

u, = V.(V u ) . (7.4) 

Let the domain of calculation be divided into an unstructured (or structured) mesh of quadrilaterals E J 

where j = 1 ,..., J . Dropping the subscript j , let the mapping from the biunit square I ■ I = I 2 = [— 1, l] 2 
onto E be denoted by 

x(^,rj) = (x(^,T)),y(^,T))). (7.5) 

This mapping can be isoparametric and the cell E can have curved edges. For simplicity, however, we 
assume that the mapping is bilinear: denoting the four corners of E by (counterclockwise) jq = jc( — 1, — 1) , 
x 2 = Jt(l,-1) , jc 3 = x(l, 1) , and x 4 = x(-l,l) , then 


x(£,ri) = (i/4) { (i - £)(i - n)x x + (i + £)(i - n)x 2 + (i + £)(i + + (i - D(i + n)x 4 }• 

See Fig. 7.1. Under this mapping, (7.3) is transformed to 


where 


J = 


+ ££ + £l. 0 

dt; dr] 


x^y n -x n y ? , u=Ju, f = y tl f-x t] g, and g= - y^f + x^g . 


(7.6) 

(7.7) 

(7.8) 


Flere, J represents the volume element. As for / (and similarly for g ), along a vertical edge, say, x x x A , 
the covariant base vector is (x ,y n ) . The vector (y -x ) represents the product of the unit normal with 

the edge length of .v, x 4 ; therefore, / = y f - x :/ g represents the normal flux across this edge (with the 

edge length accounted for). That is, the above differential formulation can be considered to be equivalent to 
the integral formulation. 

Similar to the discussion before (7.3), the flux vector / consists of two parts: f a and f . The 


advective part f a is dealt with as in (Huynh 2007). As for the diffusive part, once we obtain the solution 
for the diffusion equation (7.4), with appropriate and trivial modifications, the solution for (7.3) follows. 
Therefore, we focus only on the scalar equation (7.4) with 

f = u x ,g = u y , (7.9) 

and 

7 = y n f - X i,g = jyu - x n u y ’ and S = - Ti '/ + * *g = -y% 1 < X + ^ u y • (7. io) 

Let the solution points on I = [-1,1] be 2/t , k = 1 . For simplicity (and not necessity), the solution 

points in the i] direction is assumed to be the same as those in the £ direction: r] t = , / = 1 . The 

K- K solution points on the biunit square are (£ k ,rji ), k,l = 1,..., K . The solution points on E are 


x kj=( x k,uyk,i) = x (Zk, r h)- (7.11) 

In addition, let the solution points along the edges, which are the small squares in Fig. 7.1, be called flux 
points. For the methods discussed here, flux points always lie on cell edges. On each cell E , let the 
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solution u be approximated by K - K pieces of data u kl at the solution points x k l , k, l = l,„.,K . 
Again, boundary conditions are assumed to be periodic and are omitted. 



3 


1 ) 



(b) 


Figure 7.1. (a) Solution points (circular dots), which are Gauss points here, and flux points (squares along 
the edges) for K = 3 on the biunit square 1 2 . At each flux point, we need the common value and common 
gradient, (b) The corresponding quantities for cell E in the physical domain. 


At time level n , suppose the data u k l are known for all cells and all k, and /. For each cell, we wish to 
calculate du k ft)/dt at time t = t n (denoted by du k , / dt ) in terms of the data. That is, we wish to 
calculate u xx and u yy at the solution points. Then, we march in time by, say, a Runge-Kutta method. 

The first task is to approximate the solution by a polynomial. To this end, let <f> k denote the ID basis 

function on / given in (2.6). Let the 2D basis functions on 1 2 be defined via tensor product: for each 
(k,l) , k, l = , 

k,i = Aj(An) = ■ ( 7 . 12 ) 

Thus, (f> k 1 takes on the value 1 at (^ k ,rji) and 0 at all other K 2 - 1 solution points as shown in Fig. 7.2. 
Using the local coordinates (£, 7 ) for the cell E , the polynomial interpolating u k ] at is given by 

K K 

u E (.Av)= = Yj U k,iA(^)A( 7 l) • (7.13) 

k,l=\ k,l= 1 

Note that for each fixed 77 , u E {^^rf) is a polynomial of degree K- 1 in and vice versa. That is, u E is 
in Pk-i,k-\ ~ Pk-\ ®Pk-\ > e -S-’ for K = 2 , u E is bilinear; for K = 4 , bicubic. 


</> 


I 


(b) 

Figure 7.2. Basis functions for K = 3 with Gauss points as solution points, (a) The functions (j> 2 and (j> 2 for 
the ID case, (b) The function (j) 2 3 = f or ^ ie 2D case. 




For each fixed r/j (Fig. 7.1(a)), the interpolation along the £ -direction reduces to a ID interpolation of 
the data u k , where k varies. (A similar statement holds for each fixed f k .) On E , we can evaluate 
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(u E )e and {u E ) n (or the uncorrected gradient) at the solution points by the ID procedure. We can also 
calculate the values of n E at the flux points on the cell edge of E in a ID manner. 

Common interface values and corrected / and g at the solution points. On an arbitrary edge, say, 
in the local description, a vertical edge (similarly for a horizontal edge), let E and F denote the two 
adjacent cells as shown in Fig. 7.3(a). To define the common value and common derivative, let Xj be a 

flux point on this edge, e.g., x f = x E (l, r/,) = x F (-l, rj,) . At this point, with u E (Xf) = u E (\,t) I ) and 
u F (Xj ) = z/ f (- 1,77 ; ) , let the centered common value be defined by 

U com = u com ( Xf ) — (1 /2)[u E (x f ) + U F (x f )\. (7.14) 








E 

F 

x f 







Figure 7.3. (a) Employing u^ B (^) (corrected for the right boundary’) and u F B (^) (corrected for the left 
boundary) to calculate the common gradient at x f results in a formula of two-cell stencil, (b) Employing 
u E and u E to calculate the common gradient results in a formula of eight-cell stencil. 


The one-sided common value can be defined in a manner similar to (3.1) provided that a selection of 
left and right sides is assigned to the interface. Next, on a general cell, say, E , with the common values at 
all flux points defined as above, we can obtain the corrected £ -derivative as in the ID case of (3.10): for 
each fixed i ), , 

(«f)f(M = iu E )^n,) + [« com (-i, n,)-u E (- 1, >7/)Kg LB )'(£) 

(7.15) 

+ [ M com (l,/ 7/ )-%(l,/7,)](g RB )'(^) 

The corrected // -derivative is calculated in the same manner. At the solution points, the values of the 
corrected £-and rj -derivatives, namely u? = (u E )A% k ,Th) and u n =(u E ) tl {^ k ,rij) can easily be evaluated. 
The quantity Vz/ follows by the chain rule: 

u x = 4 + u n n x an d u y = u d y + u rj ')y ■ (7.16) 

Then, the corrected / = y l; u x - x t .u , and g = —y^ u x + x E u Y can be obtained at all solution points of E . 

Common gradients. We now define the common gradient or, more precisely, in the local description, 
/ com on the vertical edges and g com on the horizontal ones. Again, on an arbitrary edge, let E and 
F denote the two adjacent cells as shown in Fig. 7.3(a). At a flux point x , on this edge, similar to the ID 

case, if we employ u E and it ( j, to calculate the common gradient, the resulting formula generally has a 
large stencil involving eight cells as shown in Fig. 7.3(b). 

A formula for / com that has a compact stencil involving only the two adjacent cells E and F can be 
obtained as follows. Without loss of generality, assume again the edge is vertical, i.e., for some fixed /// , 
Xf = x E (\, rjj) = * F (-1, r/i) (see Fig. 7.3(a)). We can define u E (^) (corrected for the right boundary of 
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E ) and (corrected for the left boundary of F ) in a manner similar to (3.7) and (3.8). For the cell 

E , the corrected £ -derivative at x, is 

U£,E=("E)f B m = (^) f (l) + [» com -^(l)][(g RB )'(l)] (7.17a) 

and, for the cell F , 

u s , F =(u F )f(- 1) - (^) f (-l)+[M com -^(-l)][(g LB )' (-!)]• (7.17b) 

We now discuss the evaluation of u lf at x f . Since the common values u com at the flux points along 
the common edge are readily available, they can be employed to evaluate u n . Alternatively, we can 
calculate u with no correction by first obtaining the values of u E at the flux points along the common 
edge (and then u F ); however, this procedure is slightly costlier. Note that for the case of a uniform 
Cartesian mesh (used in Fourier analysis), z/ (; at x f is not needed. 

Forced E , with M;=(m £ ); ( 1) via (7.17a) and u n discussed above, we can evaluate (Vz/) £ at x f by 
the chain rule. The normal flux along the common edge namely ( f) E can then be estimated by the formula 
/ = y n u x - x^Uy . Similarly, in the case of cell F , with =(u F )®(- 1) via (7.17b) and u n discussed 

above, we can evaluate (Vu) F at x* by the chain rule, and then evaluate (f) k . To define / c< " n , we 

could employ the weighted average (4.3) provided a selection of left and right sides is assigned to the 
interface. For simplicity, the centered formula is shown here: 

7 co “ = (l/2)[(7) £ +(7) f ]. (7.18) 

The above 2D extension can easily be applied to the schemes of type I - continuous . More precisely, at 
a flux point x f , the quantity z/ com is considered to be an unknown, and it is calculated by requiring that 
( f) E = ( f) F . Flere, assuming that the common edge is of vertical type for cell E , the calculation of u v 
with no correction by first obtaining the values of u E at the flux points simplifies the scheme (see the 
discussion after (7.17)). 

The 2D extension for scheme I - recovery or (4.8) to the case of a uniform mesh of squares is simple; 
however, that for a general quadrilateral mesh is considerably more involved. 

Algorithm. The algorithm for the 2D diffusion equation is summarized below. Let the data u k t at the 

solution points x k t be given for all cells. 

□ On each edge, compute u com and / com at the flux points as follows. Assume that edge e is of 

vertical type (horizontal type is similar). Let E and F be the two cells sharing this edge. First, 
obtain the values of u E and u F at the flux points on e. Next, at each flux point Xj of e, 

evaluate and store u mm by the centered formula (7.14) or the one-sided formula (3.1). Using this 

z/ com , calculate the corrected (at one boundary only) £ -derivative u FE (xj) and u FF (Xf) by 

(7.17a) and (7.17b). Compute u„(x/) using the common values w com at the flux points along e. 

Using u k E {Xf) and u n (Xf ), evaluate (Vz/) £ at Xj by the chain rule. Again at Xj , estimate 

(f)E via / = y„u x — x„Uy . Evaluate (Vz/) /; and (f ) F in a similar manner. Calculate and store 

/ com by (7.18). Alternatively, u corn and / com can be obtained via the continuity requirement 

(4.7). Store the values w com and y com at the flux points of all edges. 

□ For each cell, calculate u, at the solution points as follows. Using w com on the four edges of the 
cell, compute (u E )^ and (z<^) 7 at the solution points via (7.15). Again at the solution points, 
calculate Vzz by the chain rule, and then / = y u x -x„u v and, similarly, g = -y^u x + x^u y . 
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Next, using / com at the flux points on the two edges of vertical type, obtain the corrected 
derivative (f)c at the solution points. Similarly, using g corn at the flux points on the two edges 
of horizontal type, obtain the corrected derivative ( g) lf . Finally, calculate u, at the solution points 
via (7.7). This completes the algorithm. 


Fourier analysis in 2D. We now carry out the Fourier (Von Neumann) accuracy and stability analysis 
for the 2D case. On the domain (- 00 , 00 ) • (- 00 , 00 ) , consider the diffusion equation 

u t = u xx +u yy . (7.19) 

The cells are the squares E t -=[i — 1/2,/ + 1/2] • [j -1/2, j + 1/2] centered at (/,/) . Denote by u Lj the 
column vector of K 2 components 

u i,j ~ ( u i,j , 1,1’ u i,j, 1,2’ •••’ u i,j, 1,JC> u i,j, 2,1’ •••’ lt i,j, 2,K> •••’ u i,j, K, 1> •••> u i,j, K,K ) • (7.20) 

Note that for the cell (i,j) , at a fixed ( k 0 ,l 0 ) , the evaluation of u xx is exactly the same as the ID case 
using the data (k,l 0 ) , where k varies in the cell (i,j) , ( i - 1 , /) , and ( i + 1 ,j) . A similar statement holds 
for w . Thus, the 2D solution follows immediately from the ID solution, and the result can be written as 


du , 


dl 


~ Q,0 u i,j + C.,0 u i-\,j + Q.0 U i+l,j + C 0,-1 u i,j-l + C 0,1 u i,j + 1 ■ 


(7.21) 


Flere, C 0 0 , C_ l 0 , C 10 , C 0 , , and C (l ] are A 2 - K~ matrices. Let the wave number in the x -direction 
be denoted by w x , and that in the y -direction, w (involving no partial derivative here). Let the 
imaginary number be denoted by I instead of the standard notation i to avoid confusion with the index i . 


we obtain 


Set 


j, u i+i,j j , and u iJ+1 respectively by e ' u t j , 

I W r ~t 'f'v . 

^ U iJ> e U iJ-- 

, and e IWy u i j, 

. I w x . Iw x gs . —I bVy 

, - 1 ^ 0 , 0 + e C_i o+ e C 10 + e C 0 j + 

at 

Iw v „ | 

e C 0,l )"i,j ■ 

(7.22) 

S = C 0 0 + e - /w * C _ 1>0 + C 1>0 + e~ Iw > C 0M + 

J w > r 
e ^ 0,1 - 

(7.23) 

■ K 2 matrix. The semidiscretization results in 



dll: ■ 

dl ' J 


(7.24) 


For all the schemes discussed here, S has K 1 eigenvalues. The eigenvalue that approximates the exact 
eigenvalue — there is one and only one such value — is called the principal eigenvalue and denoted by 
S(w x ,w y ) : 


S(w x ,w y ) « - w x 2 - w v 2 . (7.25) 

All others are spurious eigenvalues. All eigenvalues must lie in the left half of the complex plane for the 
semidiscretization to be stable. 

It turns out that for all methods discussed, via Fourier analysis, the 2D and corresponding ID schemes 
have identical order of accuracy. Concerning stability limits, the 2D method has a limit of 1/2 that of the 
ID version for nearly all schemes and, for a few exceptions, approximately 1/2. In other words, the largest 
magnitude of all eigenvalues for a 2D scheme is essentially twice that of the corresponding ID method. 
This reduction by a factor of 2 in the stability limit for the 2D diffusion equation here is consistent with the 

reduction by a factor of for the 2D advection equation in (Huynh2007): each derivative corresponds to 
a factor of J2 . 

The minimum eigenvalues, orders of accuracy, and errors of 16 schemes for the 2D case with K = 4 
are shown in Table 7.1. In the calculations of order of accuracy, the coarse mesh corresponds to w x ~ n ! 8 
and w v = nt 10 , and the fine mesh, w x =nl 16 and h’ v =^/ 20. Note that all schemes are of order at least 
6 or 2(K -1) , and the orders of accuracy are identical to those of the corresponding ID schemes in Table 
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6.3. Again, schemes 11 and 15 are of order 8, and scheme 14 (recovery), order 10. For scheme 6, or 
I - centered - g Ga /SP - g Ga , the error crosses zero near w x =n ! 8 and w v = nl 10. Therefore, to calculate 

the order of accuracy, we set w x =nl 16 and w v =;r/20 for the coarse mesh and w x =n!l>2 and 
w y =nl 40 for the fine mesh. The errors are respectively -4.64- 10~ 12 and -2.32- 10~ 14 , and the 
scheme is of order 6. The recovery scheme (I-recovery/SP-g DG or 14) is more complicated than the 
others, but it has the advantages of highest accuracy and largest stability limit. 


Schemes in 2D 
K = 4 

Min. 

Eigval. 

Ord. 

Ace. 

Coarse mesh error, 
w x = Jt/% 
w y = zr/10 

Fine mesh error, 
W x =7T / 16 
w y = k ! 20 

1. I- centered -g Le /SP-g DG 

-374. 

6 

-6.22- 10~ 9 

-2.53- 10" 11 

2. I- centered -g DG /SP-g DG 

-340. 

6 

-5.94- 10' 9 

-2.50- 10~ n 

3.1- centered - g Ga /SP - g DG 

-340. 

6 

-4.51- 10~ 9 

-2.33- 10 _n 

4. I - centered - g LumPi Lo /SP - g DG * 

-340. 

6 

5.25- 10~ 9 

2.05- 10~ u 

5. I- centered -g DG /SP-g Ga 

-244. 

6 

-1.18- 10~ 9 

-5.93- 10 -12 

6. I- centered -g Ga /SP-g Ga 

-196. 

6 

see explaination 

-4.64- 10 -12 

7. I - centered - g Lump , Lo /SP - g Ga * 

-196. 

6 

7.76- 10~ 9 

3.03- 10~ n 

8.1- centered - g DG /SP - g Lump Lo 

-244. 

6 

2.54- 10~ 9 

8.89- 10- 12 

9. I - centered -g Ga /SP-g Lump Lo 

-177. 

6 

3.42- 10~ 9 

9.91- 10- 12 

10. I - centered - g LumP; Lo /SP - g LumPi Lo * 

-165. 

6 

9.85- 10~ 9 

3.84- 10^ n 

11. I- continuous -g Le /SP-g DG 

-367. 

8 

9.50- 10~ 13 

9.28- 10^ 16 

12. I- continuous -g DG /SP-g DG 

-244. 

6 

2.62- 10~ 9 

1.02- 10 _n 

13. I- continuous -g Le /SP-g Ga 

-340. 

6 

4.91- 10~ 9 

1.92- 10" 11 

14. I-recovery/SP-g DG 

-135. 

10 

-1.89- 10~ 14 

-4.67- 10~ 18 

15. I - one - sided - g DG /SP - g DG 

-878. 

8 

2.17- 10- 12 

2.12- 10- 15 

16. I - one - sided - g Lump , Lo /SP - g Lump , Lo 

-544. 

6 

1.75- 10~ 8 

6.83- 10" 11 


Table 7.1. Minimum eigenvalues, orders of accuracy, and errors for schemes in 2D with K = 4 . 


8. Numerical Results 

As preliminary tests, we apply the schemes discussed above to solve the steady diffusion equations (or 
two-point boundary value problem). Here, the problem is solved not by time-marching, but by inversion of 
the implied coefficient matrix of u j k , j = \,...,J and k = \,...,K. Unless otherwise stated, the solution 

points are the Gauss points and, in the graphs, the crosses represent the exact solutions, while the dots, the 
numerical ones. 

The first test is the Poisson equation 

“«= 2 ( 8 . 1 ) 


on the domain [-1,1] with boundary conditions 


«(-!) = «( 1 ) = 1 . 


( 8 . 2 ) 
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The exact solution is u exact (x) = x . We employ piecewise linear methods to solve this problem. On each 
cell Ej of center Xj and length hj , the two (Gauss) solution points are denoted by x f j and Xj 2 . The 

projection of x 2 onto the space of linear functions on Ej is 

lj(x) = x 2 + hj 2 / 12 + 2 Xj ( x - Xj) . 

That is, lj provides the best linear approximation to x 2 on Ej in the least-squares sense. Note that l j is 

identical to the line defined by the values of x 2 at the two Gauss points as shown in Fig. 8.1. (In general, if 
p is of degree K , then the projection of p onto P K 1 ( /T ,) is determined by the values of p at the K 
Gauss points.) It appears desirable that the linear solution recovers f : at the Gauss points, 

«/, l = (*/, i ) 2 and u j, 2 = (*/, 2 ) 2 • 

However, as shown in Fig. 8.1, only two schemes yield this solution: I -continuous -gLe/SP-g DG and 
I - recovery/SP - g GG • All other schemes, e.g., BR2 and CDG do not. Concerning Fig. 8.1(c) for the CDG 
scheme, note that for each cell except the last one, the value of the solution polynomial at the right 
interface — thus, the common value — is exact. That is, at x - -1 , i.e., the left boundary, zqj™ = u ( /2 = 1 ; at 

each interior interface, Uj+U 2 = u j+\n = l (/0) • At the right boundary, however, the choice is switched: 
u 3 +i /2 = u 3 + 1/2 ~ 1 • This switch causes a loss in accuracy: loosely put, the cancellation of errors no longer 
holds. Also note that for all solutions of Fig. 8.1, the common interface values i/ y c °“ 2 are exact. Concerning 
the cell average values, however, the BR2 and CDG schemes do not recover the exact solutions. 




and I - recovery/SP - g DG I - centered - gD G /SP - g d G I - one sided - g ixi 'SP - gDG 

Figure 8.1. Solutions for u xx =2 on [—1,1] with boundary’ conditions w(— 1) = zz(l) = 1 using three cells. 
Only the schemes I - continuous - g Le /SP - g DG and I - recovery/SP - g DG recover the projection of the 
exact solution « exact (x) = x 2 onto the space of piecewise linear functions (the lines through the crosses). 


The second test is the Poisson equation employed in (Van Leer and Nomura 2005) 

= s(x) = -4;r 2 sin(2;Dc) (8.3) 

on the domain [0,1] with boundary conditions 

m a .( 0) = 2;r-l and z/(l) = 0. (8.4) 

The exact solution is w exact (*) = sin(2^x) + 1 — x . The projection of s(x) onto P K ~ x (Ej) is denoted by 
s j(x) . At the solution points, the source terms are set as, for k = 1 ,...,K , s - k = Sj(xj k ) . 

Figure 8.2 shows the solutions with K = 2 using three cells: in the order of decreasing accuracy, (a) 
I - continuous - g Le /SP - g DG , (b) I - centered - g DG /SP - g DG or BR2, and (c) I - one - sided - g| )G /SP - g DG 
or CDG. The solution by I - recovery/SP - g DG is nearly identical to that by I - continuous - g Le /SP - g DG 
and is omitted. The smooth green curves are the exact solutions and the dots, the approximate solutions. 
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The brown curves are the reconstruction polynomials { u C j } ; for all solutions, [i/j } are more accurate than 
the piecewise linear solutions. Concerning Fig. 8.2(c) or the CDG solution, at each interior interface, 
U j+V2 = u j+i/2 = u j( 1 ) j but at the right boundary, the choice is switched: «3°“ 2 = *<3+1/2 = 0 . Note that this 

solution is asymmetric. Also note that «y+“ 2 t°r all solutions are again exact. 


1.5 

1 

0.5 


-0.5 

(a) I - continuous - g^/SP - g DG (b) Scheme BR2 (c) Scheme CDG 

Figure 8.2. Solutions for u xx = — 4/r 2 sin(2/zx) on [0,1] with boundary conditions i/ v (— 1) = 2/r - 1 and 
i/(l) = 0 using three cells. 




Scheme 

(K=4) 

Order of 
accuracy 

Max. Error; 
4 cells 

Error 

ratio 

Max. Error; 
8 cells 

Error 

ratio 

Max. Error; 
16 cells 

I - recovery/SP - g DG * 


.215- 1(T 5 


.712- 10~ 8 


.280- 10- 10 


8 (or 2K) 


302. 


236. 


I - recovery/SP - g DG 


.769- 1(T 4 


.302- 10- 5 


.100- 10^ 6 


5 (orW + l) 


26. 


30. 


I - continuous - g Le /SP - g DG 


.433- 1(T 4 


.170- 10~ 5 


.562- 10~ 7 


5 (or K + 1) 


25. 


30. 


I- centered -g DG /SP-g DG 


.115- 1(T 2 


.106- 10- 3 


.805- 10~ 5 


4 (or K) 


11. 


13. 


I - one sided - g DG /SP - g DG 


Lh 

O 

1 

to 


.877- 10^ 4 


.574- 10~ 5 


4 (or K) 


13. 


15. 


I - centered - g Gauss /SP - g L um P , lo 


SO 

O 

1 

to 


.151- 10~ 3 


.912- 10~ 5 


4 (or K) 


12. 


16. 



Table 8.1. Maximum errors and corresponding orders of accuracy. For the line I - recovery/SP - g DG * the 
exact values at the solution points in the error calculations are the values of the projection of the exact 
solution onto P K ~ l (E •) ; for all other cases, the exact values (with no projection ) are employed. 


Table 8.1 shows the maximum errors at the solution points and orders of accuracy of various schemes 
with K- 4 for (boundary value) problem (8.3)-(8.4). As for the recovery scheme, we consider two types 
of errors: in the first type, the exact values at the solution points are the values of the projection of the exact 
solution onto P KX (E ■) ; in the second, the standard exact values (with no projection) are employed. Note 

that with the first type of errors, the recovery scheme is of order 2 K ; with the second, order K + 1 . For all 
other schemes, the two types of errors are comparable; therefore, only errors of the second type, which are 
more straightforward, are listed. These results on order of accuracy were tested for K from 2 to 7. With 
K = 2, Fourier analyses showed that most schemes are of order 2(/f - 1) = 0 ; for the steady state case 
considered here, however, these schemes are of order 2. Although not listed, all interface quantities as well 
as cell average values are exact (up to machine errors). 
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The third test is a 2D problem: on the domain [0,1] • [0,1] , with a = .2 , consider the Poisson equation 

u xx +«yy =s{x,y) = - sin(3;r x) sin(3;r y) . (8.7) 

a 

The Dirichlet boundary conditions are 

u{0,y) = w(l, y) = u(x, 0) = u(x,\) = 0 . (8.8) 

The exact solution is 

“exact (x,y) = - 1 2 sin(3^x) sin(3^y) . (8.9) 

2(3^a) 

The domain [0,1]- [0,1] is divided into J ■ J uniform squares. On each cell (square), we employ K ■ K 
solution points, which are the Gauss points. For simplicity, at the solution points, the source term is 
evaluated with no projection: for i,j = and k,l = 1 , 

s i,j,k,l = s ( x i,j,k,l) • ( 8 . 10 ) 

Figure 8.3 shows (a) the exact solution and (b) the solution with 4- 4 cells and K = 4 using 
I - centered - g DG /SP - g DG or BR2 scheme. The solutions by other schemes are very similar. 




(a) Exact solution 


(b) Solution by BR2 scheme 


Figure 8.3. Solutions for problem (8.7), (8.8) on [0,1]- [0,1]: (a) exact and (b) with 4- 4 cells and K = 4 


using I - centered - g DG /SP - g DG . 


Scheme 
(K = 4) 

Order of 
accuracy 

Max. Error; 
2 cells 

Error 

ratio 

Max. Error; 
4 cells 

Error 

ratio 

Max. Error; 
8 cells 

I - recovery/SP - g DG 


.358- 10~ 2 


.132- 10~ 3 


.378- 10~ 5 


5 


27. 


35. 


I - continuous - g Le /SP - g DG 


.237- 10~ 2 


.720- 10^ 


.214- 10- 5 


5 


33. 


34. 


I- centered -g DG /SP-g DG 


.117- lO^ 1 


.851- 10~ 3 


.116- 10~ 3 


» 4 


14 


7 


I - one sided - g DG /SP - g DG 


.117- 10~ ! 


OO 

O 

O 

1 

TJ 


.127- 10- 3 


» 4 


7. 


14. 


I - centered - g Gauss /SP - g LumPi Lo 


.450- lO^ 1 


.147- 10~ 2 


.134- 10~ 3 


« 4 


30 


11 



Table 8.1. Maximum errors of values at solution points and corresponding orders of accuracy. 
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Table 8.2 shows the maximum errors at the solution points and orders of accuracy of various schemes 
with K - 4 for the 2D problem (8.7)-(8.8). Note that the orders of accuracy are essentially the same as 
those of Table 8.1. 


9. Conclusions and discussion 

In summary, a new approach to high-order accuracy for diffusion problems with the advantages of 
simplicity and economy was introduced. The approach evaluates the derivative for the discontinuous 
solution polynomial by first obtaining the derivatives at the solution points in a straightforward manner. 
These derivatives are then corrected, in each cell, by quantities proportional to the jumps at the two cell 
interfaces. The key idea is to solve the equations in differential form using a reconstruction technique 
where the discontinuous solution polynomials are approximated by piecewise polynomials that are 
continuous across the interfaces and of one degree higher than that of the solution polynomials. The 
approach resulted in several new methods including a simplified version of the DG scheme. It also led to 
new choices of common value and common gradient. Fourier stability and accuracy analyses were carried 
out. For the two-point boundary value problem, it was shown that nearly all of these reconstruction 
schemes (which include a majority of popular DG methods) yield exact common interface quantities as 
well as exact cell average solutions. 

The approach can be extended to the case of the Navier-Stokes equations on a 2D quadrilateral mesh 
via tensor products. The extension to the case of a triangular mesh appears to be more complex and remains 
to be explored since the concept of tensor product is no longer available. 
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Appendix. Boundary Value problem (or Steady State Case) 

In this appendix, we show that for the ID boundary value problem, assuming that the solution exists, 
then the common interface values and derivatives and the cell average values of the solution are exact for 
essentially all schemes discussed. 

More precisely, on the domain [a, b ] , consider the Poisson equation 


U xx=* 

with the Neumann boundary condition on the left and Dirichlet boundary condition on the right, 

u x (a) = a and u(b) = /3 


(A.1) 

(A.2) 


where a and /? are given constants and, for simplicity, s = s(;c) is assumed to be piecewise continuous. 
Let "exact be the exact solution and v exact = (w exact ) z . Then, for y in [a, b \ , 


Vexact (y) = a+ [' S(z)dz 

Ja 

(A.3) 

and, by integrating from right to left. 


"exact (-*■) — P + j* ^ v exact (y ) dy . 

(A.4) 

Next, consider the following three sets of boundary conditions: 


II 

$ 

II 

(A. 5 a) 

u(a) = a, u x (a) = a, and 

(A. 5b) 

s 

o- 

II 

II 

(A. 5 c) 


We now show that each of them can be recast in the form (A.2). For the case (A. 5a), by (A. 4) and (A.3), 


"exact {x) = P + a(x-b) + f [ s{z)dzdy. 

Jb Ja 

Therefore, 

"exact {a) = p-a(b-a)~ f [ s(z)dzdy . 

Ja Ja 


Or 


a - P - a{b- a)- \ s(z)dzdy . 

*a Ja 


If a and ft are given, we can evaluate a in terms of a and fS by the above equation. Therefore, (A. 5a) 
can be cast in the form (A.2). The other two cases are similar and are omitted. 

On domain 7 = [-1,1]. First, we consider the case of only one cell on the domain [a, A] = [-1,1] . 

Assume that K> 1 . Denote by S' the projection of s onto P K l (I ) , the space of polynomials of degree 
K - 1 or less on I . Using P k , the Legendre polynomial of degree k , set 

tr-i 

? = P r A_iW = J^r k p k 

k = 0 
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where, for k = 0 ,..., K - 1 , y k is defined by 

Yk( P kA) = (S’ P k)- (A.6) 

The exact solutions for the problem (A. 1) and (A. 2) on the domain [-1,1] are 

VactO) = «+ (A. 7) 

and 

^exact (-*0 — P ^exact (T) ^ • (A.8) 

The exact solutions for the same problem with s replaced by s are 

v(y) = d+ (A.9) 

and 

u(x) = /? + v(y)dy . (A.10) 

Here, since 7 is of degree K - 1 , v is of degree K , and u , degree K + 1 . 

We claim that 

J ■ equivalently, v(l) = v exact (l) . (A.lla,b) 

In addition, if K> 2 , then 

\j({)d{ = j_v exact (t;)A ; equivalently, z7(-l) = w exact (-l) . (A.12a,b) 

And, if K > 3 , then 

\_u(Z)d% = {^exact (^)d^ ■ (A. 13) 

To prove this claim, note first that the projection of 7 and s onto P K 1 ( / ) are equal; that is, for any 
polynomial p of degree K - 1 or less, 

(s,p) = (s, p ) . (A. 14) 

If K > 1 , then P K ^(I) contains p = 1 . The above for p = 1 yields 

[s{g)d$=\[s{£)d$, (A. 15) 

which is (A. 11a). By (A.7) and (A.9), equation (A. 11a) is equivalent to equation (A. lib). 

Next, if K > 2 , then P K ~ l (I) contains p(, g) = ^ . Equation (A. 14) for this p yields 

jy (£)£</£= jy (#)£«/£. (A. 16) 

Note that by (A.9), v '(£) = s(£) . Applying integration by parts to the left hand side above, we obtain 

L/^di= J_jV '4^ = v(l) + v(-l)- \ vd£. 

As for the right hand side of (A. 16), 

j* j ~r d ~ — J ^ v cxac t 7 dc — V exac t(l) "t" H-xact ( — ^) — J | A'xact d% . 

Since v(-l) = v exact (-l) = cc and, by (A.l 1), v(l) = v exact (l) , the above three equations imply 

J_i vd£= J ( v cxact d % , 

which is (A.12a). By (A.8) and (A.10), equation (A.12a) is equivalent to equation (A.12b). 


(A 17) 
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Finally, if K > 3 , then P K l (I ) contains p = . Equation (A. 14) for this p yields 

Jy(£)£ 2 </£ = Jy(£)£ 2 </£. (A. 18) 

Again, by (A.9), v '(£) = s(£) . Applying integration by parts to the left hand side above, we have 

J^s^ 2 = J 1 v '^ 2 dt; = v(l) - v(— 1) - 2 ^ v %cU; . (A.19) 

In a similar manner, by (A. 10), u'=v. Applying integration by parts to the last integral above, we obtain 

J v^= J «'^ = «(l) + af(-l)- ^udl;. (A. 20 ) 

Substituting the above into (A.19), i.e., after integrating by parts twice, the result is 

| s£ 2 dl; = v(l) - v(— 1) - 2w(l) - 2m(— 1) +2 J ud^ . 

A similar argument for the right hand side of (A. 18) leads to 

fl 2 fl 

— ^exact(^) — ^exact ( — 0 — ^^exact(^) — ^^exact ( — ^) + ^ J ^exact • 

Recall that that v(-l) = v exact (-l) = a and, by (A. 11), v(l) = v exact (l) ; in addition, w(l) = M exact (l) = /? 
and, by (A. 12), w(-l) = » ex act( - 0 • Therefore, the above two equations and (A.18) imply 

j U(£)d%= J^exactC 1 ?)^- 

This completes the proof. 

Note that an argument along the same line as the above for the recovery scheme was discussed in (Van 
Leer et al. 2007). 

Exactness of common interface quantities and cell average values regardless of the scheme. Let the 

domain \a,h\ be partitioned into (possibly irregular) cells Ej of center x } and length hj, j = 1, ...,./ . 

Assume that K> 1. Let x j k , k = l,...,K, be any set of solution points. On each cell E-, denote by 
Sj=Sj(x) the polynomial of degree K - 1 that is the projection of s onto P k ~ x {Ej)\ 

*j = FV-i(s) • (A.21) 

For the s' values at the solution points, we employ, with k = 1 , 

s j,k= s j( x j,k) ■ (A. 22) 

The problem takes the form: find u , k such that for all j and k 

( u xx)j,k =S j,k 

with the following boundary conditions: 

( u x)\/2 ~ ( w x)l/2 ~ a ar *d lt J+U2 ~ U J+l/2 ~ P ' 

At each interface j + 1/2 , let 2 be the common value and (M v )y°“ 2 the common derivative. They can 
be calculated from u. k and w + lk , k = l,...,K. 

The following claim is in order. Assuming that the solution exists, then, regardless of the method, (a) 
the common derivatives are exact; (b) if K > 2 and the correction function g for the solution points has a 
zero average on [-1,1], then the common interface values are exact. In addition, (c) if K > 3 and the 

correction function g for the solution points is orthogonal to P 1 , then the cell average values of the 
solution are also exact. 
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Concerning statement (b) in the above claim, since g is orthogonal to P K 3 , g has a zero average if 
K> 3 . Moreover, if g = g DG , then, since g DG is orthogonal to P K ~ 2 , it has a zero average if K> 2 . 
Concerning statement (c), since g is orthogonal to P K ~ 3 , g is orthogonal to P 1 if K > 4 . Moreover, if 

g = g DG , then, since g DG is orthogonal to P K ~ 2 , g is orthogonal to P 1 if K > 3 . 

To prove the claim, we employ the reconstruction functions. Recall that the reconstruction scheme has 
the following properties: 

u C j (x) is of degree K and takes on the value u “” 2 at x y _ 1/2 and u ™” 2 at x j+ 1 / 2 1 
{vy(x) = (u C j ) x (x)} is a discontinuous piecewise polynomial of degree K - 1 ; 

(x) is of degree K and takes on the value Vy°“ 2 at x y _ 1/2 and v““ 2 at Xy +1/2 ; 

\{v c : ) x } is a discontinuous piecewise polynomial of degree K - 1 . 

If Uj k are the solution of the boundary value problem, then, in each cell, the corresponding (v ( / ) x , 
which is of degree K - 1 , satisfies, for k = 1,..., K , 

(v C j) x (x jk ) = s hk . (A.23) 


Since sAx) is the (unique) polynomial of degree K - 1 interpolating s jK . the above implies 

(vj ) x (x) = s'y(x) . Therefore, recalling that v“” 2 = (t< x )y°” 2 , 


Consequently, 


C / \ / \ com 

Vj (x) = (u x )j_ 1/2 + 



(Mx) 


com 
y+1/ 2 


= («*) 


com 

y—1/2 


+ 



Concerning the above integral, since Sj = P%_ 1 (5') , we have 


(A.24) 

(A.25) 


\ ~ s Ax)dx = f J+l 2 s(x)dx . 

•bj-1/2 J - x j - 1/2 


(A.26) 


We now march from / = 1 to j = J ■ For j = 1, at the left interface of C, . namely, at x l/2 =a, the 
quantity (m x )“” = « is known and is exact. We can calculate (m x )“™ by applying (A.25) with j = 1 . 
Equations (A.26) and (A.3) then imply that (w x )™“ is exact. Continuing for j = 2, ..., J , respectively, 
we conclude that all common interface derivatives («*)>+ 1 / 2 are exact regardless of the method: for all j , 


(»,)“”/2 = (Vexact) 2+ l/2- (A. 27) 

Note that in the proof, the key ingredient is (v C j ) x = s f ; we do not need v ■ . 

To show part (b), i.e., the exactness of the common interface values m“” 2 , we integrate from right to 
left and make use of the same result in the local coordinate, namely, (A. 12). Note first that by 
(A.24), Vy (x) is independent of the method chosen; v y (x) , however, is method dependent. Therefore, we 

wish to obtain w co ” 2 - w“” 2 by integrating v C j{x) instead of Vj{x) . To this end, recall that by definition, 
Vj{x) = (m^ ) x (x) . In addition, at all interfaces, u c j(Xj +l / 2 ) = m “™ 2 . Thus, for any x on E f , 

u C j (x) = u + { Vj(y)dy . 

Jx j+ 1/2 

With x = Xj_ 1/2 , 

com com C X J ~ 1/2 / \ j . . 

u j- 1/2 - Wy+1/2 = I • (A.28) 

* ,Jc y+l/2 
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In the local coordinate, Vj(^)= vT^) + c l g- LB (£) + c 2 g RB (£) for some constants c x and c 2 . Now, 
suppose the correction functions for the solution points satisfy, for g = g LB and g = g RB , 


{ , | g(^)^ = 0. 


(A.29) 


Then, 


j ' l 2v j(y)dy = J ' V1 v C j(y)dy . 

Jx j+ 1/2 Jx j+ 1/2 

As a consequence, if (A.29) holds, then, by (A.28), 

com com C x j-V2 q . , , 

u j-\n ~ u j+\n ~ I Vj(y)dy. 

"X fa-1 n 


(A. 30) 
(A.31) 


To show that these interface values are exact, we employ (A. 12). Starting with the last cell Ej , the 
boundary quantity 2 = [i is exact. The polynomial \’j of degree K plays the role of v in (A. 9). 
Therefore, (A.12) implies that wj°” 2 is exact. By using this argument successively for j = J-1,.,.,1, i.e., 
marching from right to left, we conclude that the interface values m c °“ 2 are exact provided that K> 2 and 
(A.29) holds. 

Finally, assuming that K> 3, we wish to employ (A. 13) to show that the cell average values of the 
solution are exact. The problem here is that u of (A. 13) is of degree K + 1 whereas iij is of degree K - 1 . 

To bridge this gap, we first observe that since u ( j = U j +c 3 g LB +c 4 g RB for some constants c 3 and c 4 , 
and the averages of g LB and g RB are zero by our assumption, the average of w ; on Ej is the same as that 
of Uj (degree K ). Therefore, we only need to bridge the gap between ly (degree K ) and u (degree 
K + 1 ). To this end, recall that we integrate s ; , which is of degree K - 1 , to obtain v ( j , which is of degree 
K ; here, Vj plays the role of v . Next, we integrate v • and v R to obtain i/j , and it , respectively: for 
any x on Ej , 


^ (x) = u ™? /2 + f Vj(y)dy , 

Jx j+ 1/2 

and 

w(x) = m“” 2 +[ v C j (y)dy . 

Jx j+ 1/2 

C LB RB • • LB RB 

Since Vj =Vj + c x g +c 2 g , the problem reduces to the following. With p = g or p = g , set 


p (n)= ■ (A.32) 

We wish to show that 

\P(v)dri = 0. (A.33) 

If we can show the above, then u - , and u have the same averages on E - . 

In fact, we will show a more general and simpler to state result: if p is orthogonal to P l , i.e., 

| ^p{%)di; = 0 and | ^p{^)^d^ = 0 , (A.34a,b) 

then P defined by (A.32) satisfies (A.33). Indeed, by (A.34a) and definition (A.32) of P , we have 
P(l) = 0 . Next, using integration by parts and (A.34b), 

P{S)dS = \P{J;)$\\- \[p'{Z)Zdt = P{\) + P{-\)- \\p{Z)ZdZ = 0. 

This completes the proof. 


(A. 3 5) 


